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SOLAR DYNAMIC HEAT RECEIVER 
KYUNG-HWAN AHN 
ABSTRACT 



In the present investigation a thermal design of solar receiver has been 
developed for the solutions of problems involving phase-change thermal energy 
storage and natural convection loss. Two dimensional axisymmetrical solidifica- 
tion and melting of materials contained between two concentric cylinders of finite 
length has been studied for thermal energy storage analysis. For calculation of 
free convection loss inside receiver cavity, two dimensional axisymmetrical, lam- 
inar, transient free convection including radiation effects has been studied using 
integral/finite difference method. Finite difference equations are derived for the 
above analysis subject to constant or variable material properties, initial and 
boundary conditions. The validity of the analyses has been substantiated by 
comparing results of the present general method with available analytic solu- 
tions or numerical results reported in the literature. Both explicit and implicit 
schemes are tested in phase change analysis with different number of nodes rang- 
ing from 4 to 18. 


The above numerical methods have been applied to the existing solar re- 
ceiver analizing computer code as additional subroutines. The results were com- 
puted for one of the proposed Brayton cycle solar receiver models running under 
the actual environmental conditions. Effect of thermal energy storage on the 
thermal behavour of the receiver has been estimated. Due to the thermal energy 
storagy, about 65% reduction on working gas outlet temperature fluctuation has 
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been obtained, however maximum temperature of thermal energy storage con- 
tainment has been increased about 18%. 

Also, effect of natural convection inside a receiver cavity on the receiver 
heat transfer has been analyzed. The finding indicated that thermal stratification 
occurs during the sun time resulting in higher receiver temperatures at the outlet 
section of the gas tube, and lower temperatures at the inlet section of the gas tube 
when compared with the results with no natural convection. Due to heat supply 
from the air during the shade time, minimum temperature has been increased, 
while maximum temperature has been reduced due to convection loss to air. 
Consequently, cyclic temperature fluctuation has been reduced 29% for working 
gas and 16% for thermal energy storage containment. Despite of the presence of 
the natural convection, on the other hand, the time-averaged temperatures for 
receiver components were found to be similar for two cases with/without natural 
convection (maximum difference was 1.8%). 
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CHAPTER I 
INTRODUCTION 

In engineering practice, industrial processes, and laboratory work it is often 
necessary to find the temperature distribution and heat flow in bodies under 
transient conditions. Among these transient heat transfer problems, attention 
is given for cases involved in phase change of a material and natural convection 
inside an enclosure. 

1. ) Case related to Phase change of a material: heat transfer due to solidifica- 

tion or melting of materials is of considerable importance in many technical 
fields. Recent interest in the storage of thermal energy from intermittent 
sources such as the sun, as well as ongoing interest in many processes in 
metal casting, food freezing, aerodynamic heating of high speed traveling 
objects such as missiles and space craft and freezing or thawing of soil, etc. 
has highlighted the importances of phase change phenemena. 

2. ) Case related to natural convection inside an enclosure: heat transfer due 

to bouyance-driven flow of the fluid contained in a solid boundaries is 
encounted in many engineering applications, including the storage of cryo- 
genic fluids, petroleum storage vessels , air circulation in the building. In all 
these cases, the heat loss from or to the enclosure is an important problem. 
Chapter II provides extensive literature survey for the two-dimensional 
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axisymmetrical solidification and melting problem. This is followed by a detailed 
presentation of the different numerical methods used along with results and 
comparisons for sample cases. Similarly, in chapter III, numerical derivation 
of a heating process by the natural convection is given with a related literature 
survey. In chapter IV, these numerical methods are applied to thermal analyses 
of the solar receiver of the space station freedom solar dynamic power system 
by incorporating the phase change subroutine developed in Chapter II and the 
natural convection subroutine in Chapter III to the existing receiver program, 
namely HEAP code, which will discussed later. In addition, several sensitivity 
analysis will be conducted to examine the receiver thermal performance. 

1.1 Energy Storage Solar Dynamic Receiver 

1.1.1 Background A solar dynamic power system for generation of auxiliary 
electrical power in space has been recently developed using solar energy as heat 
source. A thermodynamic power cycle for such system is the closed-loop Brayton 
cycle (see Figure 1.1). The solax Brayton heat receiver is a combination of heat 
exchanger and heat storage device which transfers energy to a Brayton cycle 
working gas during both the sun and shade period of an Earth orbit. During a 
sun period, a parabolic collector is used to focus solar radiation into the heat 
receiver through a small aperture located at the collector focal point. The solar 
input is normally designed to be greater than that required to heat the working 
gas at the design outlet temperature so that the excess solar energy input during 
a sun period could be stored in phase change material of heat storage system 
and then withdrawn by the working gas during the following shade period. 

Thermal energy can be stored by two different ways: a) Sensible heat 
storage system and b) Latent heat storage system. The sensible heat storage 




Figure 1.1 Schematic diagram of the closed Brayton cycle solar 
dynamic power conversion system. 
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systems utilize the specific heat capacities of liquids such as water or solids such 
as rocks. On the other hand the latent heat storage systems utilize the latent 
heat of phase change material(PCM). The energy storage systems using phase 
change materials are particularly attractive for the following reasons: 

(1) Their ability to store or release large amount of thermal energy in the form 
of latent heat, per unit weight of storage material, resulting in smaller ther- 
mal energy storage system weight as compared with sensible heat storage 
systems; 

2. Their ability to store or release thermal energy during phase change at a 
nearly constant temperature. 

1.1.2 Heat Storage Receiver Energy Balance In general, the temperature dis- 
tribution of the heat receiver and of the receiver outlet gas will vary in time 
depending on the mass flow rate and inlet temperature of the gas, on the solar 
input rate, on the rate of energy losses through the aperture and through the 
insulation, on the lengths of the sun and shade periods for one orbit cycle. 

If the orbit cycle is such that the excess solar energy input available for 
storage during a sun period is much in excess of the energy requirements during 
the following shade period, then all the phase change material would melted 
and temperature become too high. Accordingly the system will be overheated 
and burn out. In order to prevent that situation, it is necessary to determine 
proper values for the aforementioned parameters which brings stable system 
performance. 

A design objective for the solar receiver is to minimize time variations in 
gas outlet temperature. Therfore, a heat balance must be carefully performed on 
the total receiver during an orbit. Figure 1.2 shows the different type of energy 
involved in the solar receiver during the sun and shade period of an orbit. These 
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energies axe: 

Qs > the rate of solar energy input; 

Q a , the rate of energy extracted by gas; 

Ql, the rate of energy loss; 

QsTi the energy stored during a sun period. 

An energy balance for these quantities during the sun period provides the 
energy that should be picked up by the PCM. In turn this energy would be 
drained from the PCM during the following shade period. 

In the general situation each of the above quantities is a function of time. For 
cyclic equilibrium operation, the energy stored during a sun period equals the 
energy extracted from the receiver during the following shade period; also equals 
energy absorbed or released by the PCM, that is, 


/ [<2s(t) - Qg(<) - QL(t)]dt = Qst = [ [Qg(0 + Ql(<)1 * 

J At^sun *■ ^ J At, shade ^ ^ 


(i.i) 


The solar energy input term (Qs) of Equation 1.1 is a function of the following: 
1) receiver and concentrator geometry; 2) receiver and concentrator radiative 
surface properties; 3) receiver-concentrator orientation with respect to the sun 
and with respect to each other; 4) solar constant and solar collimation angle. 

The heat loss term (Ql) in Equation 1.1 consists of loss of radiative energy 
through the aperture (Qa) a-nd energy losses through the insulated surfaces of 
the receiver (Qj)‘, that is, 


= J [<rTi(i)‘ - <rll] + „,((()} iAi 


( 1 . 2 ) 
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(b) 


Figure 1.2 Orbital receiver thermal energy balance, (a) During 
sun period, (b) During shade period. 
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where the integration is carried out over all of the receiver inside surface and T a 
represents the total radiation view factor from the node i to the aperature. The 
insulation heat loss ( q j) is dependent on local insulation properties and tempera- 
ture difference across the insulation. It should be noted that, under 1-g operation 
for the receiver, natural convection inside the receiver cavitywill contribute to 
the heat losses, and additional terms should be included in Equation 1.2 as will 
be shown later. 

Finally, the required heat transfer to the working gas (Qa) is determined 
by the Brayton cycle design operating points as follows: 

Qg(<) = m Cpg (r Go%t - T g (1.3) 

where Ta out and Ta in are the specified gas temperatures at exit and inlet to the 
receiver respectively. On the other hand the heat transfered to the gas can also 
be expressed as 


NP fl 

Qo{t) = £ / 

»=i Jo 


hirD 


M - r 0 (*)] 


dz 


(1.4) 


where the summation is to include all gas tubes and the integration is carried 
out along the length of each tube. Equations 1.3 and 1.4 can be used to yield the 
working fluid tube size and temperature as well as conductance necessary to meet 
the design requirements of the receiver. In Equation 1.1, one could substitute for 
Ql and Qg from Equations 1.2 and 1.4 respectively. Thus, the workinf fluid tube 
wall temperature is directly related to the surface temperature of the receiver 
cavity, which in turn is dependent on the solar input and the length of the sun 
and shade period. 
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Since the PCM is separating the working fluid tube wall and receiver cavity 
surface, a detailed heat transfer analysis for the PCM should be introduced to 
ensure a local heat balance at each point within the receiver. Such an analysis 
will be considered in Chapter IV. 

Before getting into the detailed anaysis for the phase change problem and 
natural convection loss, simplified receiver heat balance should be performed 
by assuming uniform temperature for the receiver materials to understand a 
relationship between the different type of energies described above. In order to 
minimize variation in the receiver and the outlet gas temperatures with time, the 
independent parameters including geometry should be chosen so as to keep the 
PCM in the two phase condition. As long as there is the phase change occuring 
within thermal energy storage container, the exterior surface of the gas tube wall 
temperature is expected to be neax the fusion temperature of the PCM. 

In order then to use Equation 1.1 for obtaining a gross estimate of the 
potential performance, one should assume that the PCM always remains in a 
two phase condition and that the temperature variations through the PCM can 
be neglected. This leads to the simplification that the surface temperatures 
appearing in Equations 1.2 and 1.4 are the melting temperature of the PCM 
(I'm) which is independent of time. Under these conditions Equation 1.1 can be 
rewritten as follows: 


( Qs -Qg - QL)&i»un = ( Qg + Q^&tihade = Q ST (1-5) 

Furthermore if one assume that the gas side conductance used is a mean value, 
independent of time, the expression for heat transfer to the gas can be written 


as 
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Qa = m Cpa [r. - To(0)] (l - e*® 1 ) 


( 1 . 6 ) 


where 

NPtrDh 

~ WC P , G 

Since the mass flow rate, inlet and outlet temperatures are specified by the Bray- 
ton cycle design points, Qg is also specified from Equation 1.3 and Equation 1.6 
yields the required gas side conductance necessary to meet the design require- 
ments when the tube surface temperature is held constant. 

The loss term appearing in Equation 1.5 consists of radiative losses through 
the aperture Qa, and of losses through the insulation Qi (Equation 1.2). Using 
the previously made assumption that the cavity temperature is the same as 
melting temperature of the PCM, the radiative losses through the aperture can 
be written as followes: 


Q A =A„(,Ti.-ali) (1.7) 

With the value of Qa specified by the Brayton cycle requirements and the value 
of Qa from Equation 1.7, Equation 1.5 can be used to solve either the required 
minimum solar energy input Qs for a given insulation loss Qi, or the maximum 
allowable insulation loss for a given solar input. Accordingly the energy storage 
capacity ( Qst )> defined as the product of total PCM mass and the heat of fusion 
of PCM, can be obtained. 

The actual sun-shade period encountered during an orbit mainly depend on 
the orbit altitude and on the angle between the orbit plane and the ecliptic. For 
a given circular orbit, the maximum shade time occurs when the intersection 
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of the orbital plane and the ecliptic plane coincides with the Earth-Sun line. 
Considering circular orbits for a given altitude there is a critical value of orbit 
inclination above which a satellite would encounter some orbits of continuous sun 
time during a year. The illumination conditions of several orbits are summarized 
in Table I. Although the receiver must be designed and insulated to operate 
under maximum shade conditions, these conditions are encountered only during 
a small percentage of orbits per year as shown in Table I. Accordingly, during 
most orbits, the required energy for a shade time is less than the maximum, and 
it is necessary to develop a method to reject some of the excess solar energy 
so that the receiver temperatures and the outlet gas temperature are kept at 
acceptable levels in the sun period. 

Recently, computer code has been developed for prediction of the thermal . 
behavior of solar receivers and was given the name HEAP, an acronym for Heat . 
Energy Analysis Program. HEAP has a basic structure that fits any general 
heat transfer problem with specific features that are “custom-made” for solar 
receivers. However, this code doesn’t include any analysis for a thermal energy 
storage, and accordingly can not be used to predict the thermal characteristics 
of the real receiver which has a thermal energy storage unit. Preliminary study 
has been conducted assuming no property variations between the liquid and solid 
phases without calculating the location of the liquid-solid interface resulted in a 
conclusion that more rigorous analysis which allows not only property variations 
between the liquid and solid phase but also interface motion is needed for the 
accurate prediction under time dependent environmental boundary condition. 
Since HEAP program is divided into many subroutines, another subroutine that 
solves phase change energy storage problem using the method described in chap- 
ter II is added into the existing program. 
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1.2 Earth-based Receiver Code 


Due to time and cost constraints it is likely that most components of the 
space station solar dynamic power system will undergo extensive earth based 
testing as opposed to orbital testing. This will expose them to gravity driven 
effects such as natural convection. This natural convection could influence the 
heat transfer in the working fluid loops. In most cases the Reynolds number of 
the flows will be sufficiently large so that the effect will be insignificant. The 
free convection effects could be expected to be important in two cases [Evans 
(1987)]. In the receiver, two regimes are considered where natural convection 
might be important: 1) In the working fluid tube flow, natural convection will be 
important only if the flow is running at low Reynolds number or undergoing phase 
change; 2) In the receiver cavity, particularly if the chamber where the receiver 
is tested is not evacuated (Figure 1.3). Under the current design conditions 
the working fluid flows are at high Reynolds number and kept in one phase 
throughout the receiver. Therefore, only natural convection in the receiver cavity 
will be considered in this analysis. 

Preliminary study has been made to examine the effect of free convection 
on the receiver utilizing an overall convection coefficient for the receiver cavity. 
One can estimate the size of the convection heat loss rate. While this overall 
effect was small, there was substantial changes in the wall temperature that might 
be enough to effect significantly the stresses that are generated in the receiver 
materials. It should be noted that the above overall convection heat transfer 
coefficient has been calculated using equations developed for related geometries 
and not from the test results on actual receivers. Therfore it is necessary that 
a detailed analysis of the convection inside the receiver be under taken to fully 
understand the convection effects in receiver performance. 
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Hitherto, thermal design of solax receiver in solar dynamic power system for 
ground test as well as in solax thermal central receiver system is commonly inves- 
tigated with approximate models where the receiver is treated as an isothermal 
box with lumped quantities of heat losses by free convection to the surround- 
ings. These approximate models, are not detailed enough to distinguish between 
receiver designs, to predict transient performance, or to accurately express the 
performance sensitivity to variations occuring in the system parameters. 

For more adequate determination of convection loss, it is necessary to de- 
velop a procedure which applies state-of-the-art numerical techniques to analyze 
the complexities of buoyance-induced flow in cavity-type receiver geometries, and 
provide finite difference solutions to the equations governing mass, momentum, 
and energy conservation for the axisymmetrical fluid flow. 

1.3 Objective and Scone 

The main objective of this thesis is to thermally analyze the performance 
of the heat receiver under different orbital and terrestrial operating conditions. 
This objective can be achieved by conducting systematic study on following: 

(1) The phase-change problem in order to find the temperature distribution 
and solid-liquid interface motion of the thermal energy storage (TES) ma- 
terial, and heat transfer rates in and out of the material in Chapter II; 

(2) The enclosure natural convection problem in order to find the temperature 
and velocity distribution of the fluid inside the enclosure and heat transfer 
coefficient in Chapter III; 

(3) The effect of (1) and (2) on the receiver thermal performances under real 
environmental conditionsin Chapter IV. And the above analysis will be 
su mm arized and resulting conclusions will be drawn in Chapter V. 
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At this point it is obvious that the natures of these transient problems are 
such that no known general solutions exist and the problems must be solved by 
using numerical techniques. In the present thesis, therefore, efficient numerical 
scheme will be developed and presented with comparison with available data and 
finally applied to the solar dynamic receiver thermal analysis. 
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CHAPTER II 

HEAT CONDUCTION ANALYSIS OF PHASE CHANGE PROCESS 

2.1 Introduction 

Heat conduction problems involving freezing or melting axe complicated 
due to the coupling of the temperature field with the rate of propagation of 
the phase boundary between the solid and liquid phases. The first published 
discussion of such problems seems to be that by Stefan in a study of the thickness 
of polar ice [1891], and for this reason the problem of freezing is frequently 
referred to as the problem of Stefan. 

The essential new feature of such problems is the existence of a moving 
surface of separation between the two phases. Associated with this, the way in 
which this surface moves has to be determined. Heat is liberated or absorbed on 
it, and the Thermal properties of the two phases on different sides of it may be 
different, so that the problem is highly non-linear and special solutions must be 
determined and cannot be superposed. 

In the past, there has been many attempts to obtain solutions to these 
problems. Only for one dimensional problems, however, fruitful results were 
obtained. This is due to the complexity of the problem whose solutions have to 
be obtained simulteneously from the following three non-linear partial differential 
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equations [Ozisik (1980)]: 



= V • (*, VT,) 


(2.1a) 


k t VT, • VF - kiVTi • VF = - pL ^ 

Ot 


(2.16) 


C Pl p-^ = V'{k l VT l ) (2.1 c) 

where Equation 2.1a and 2.1.C axe the conduction heat equations for solid and 
liquid phases respectively. While Equation 2.1b is the interface condition to be 
satisfied at the surface of separation. 

In most cases idealizations and simplifications have been made to obtain 
solutions for the above general equations. In Figure 2.1, three-dimensional solid- 
ification is illustrated where the solid and liquid phases are separated by a sharp 
interface defined by the equation 

F(x,y,z,t) =0 

and n denotes the unit vector normal to this interface at any location P on the 
interface and pointing toward the liquid region. 

A survey of the literature shows the major assumptions on which these past 
solutions are based and indicates the need for future research. Only a few exact 
analytical solutions have been found and most available solutions are obtained 
by analytical approximations and numerical methods. The principal approachs 
to a solution of the solidification processes which have been tried in the past can 
be divided into two classes: 



Figure 2.1 Solidification in three dimensions: interface moving in 
the n direction. 
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(a) Analytical method 

(b) Numerical method 

Each of these groups will be discussed separately. 

(a) Analytical Methods 

For constant material properties and for proper initial and boundary con- 
ditions formal analytical solutions can be obtained. In most cases, however, the 
results are seriously limited, on one hand by the assumptions made, and on the 
other hand by the difficulties encountered in finding proper solution methods. 

The most well-known exact solution is that of Neumann for the semi- 
infinite region x > 0 initially at a constant temperature which is greater than the 
melting-point, and with the surface x = 0 maintained at melting temperature. 
Exact solutions are also available for a number of problems on the infinite region 
in which x < 0 is initially solid at constant temperature and x > 0 is initially 
liquid at constant temperature. 

No exact solutions are available for other geometries such as the slab —a < 
x < a with its surface maintained at zero, or for the region —a < x < a initially 
liquid and |z| > a initially solid. 

One-dimensional problem for simple shapes has been treated by various 
analytical techniques, e.g., the perturbation method by Jiji (1970), the heat 
balance integral method by Goodman (1958), the variational method by Biot 
(1957), the method of moving heat sources by Lightfoot (1929), and the method 
of polynomial approximation by Megerlin . Analytical solutions are scarce for 
problems of more than one dimension. 

For problems on radial flow in cylinder or spherical coordinates, the only 
simple exact solution in cylindrical coordinates corresponds to supply or removal 
of heat by a continuous line source. An exact solution is given by Ingersoll for 
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materials extended radially to infinity and cooled by a line source with constant 
heat flux. For the region bounded internally or externally by a circular cylinder 
with constant surface temperature, only an approximate solution is available. 

London and Seban obtained solutions of quasi-steady temperature distri- 
bution to problems involving freezing processes of water. By assuming that the 
heat capacity of the solidified or melted substance is negligible relative to the 
latent heat of fusion, they deduced exact closed-form solutions. The equations 
resulting from this method describe the fusion front motion in bars, cylinders 
and spheres for different heat transfer coeffcients but are restricted to cases 
where the materials extend to infinity, where the initial temperatures are at the 
freezing point. They concluded that an analysis of the degree of approximation 
introduced by the idealization of negligible capacity effect indicates that the re- 
sulting solutions are sufficiently accurate for most engineering calculations and 
the neglect of ice thermal capacity introduces no greater error than the uncer- 
tainty of the correct magnitude of the unit thermal conductivity k. This analysis 
was further expanded for general orthogonal curvilinear coordinates and validity 
of neglecting the sensible heat effects for low Stefan numbers was confirmed by 
Bledjian (1980) who compared his results with Magerlin (1968) and S hamsunder 
and sparrow (1975). 

Most of the analytical investigations presented above neglect latent heat 
effects and the thermal capacity of the solid and doesn’t include finite heat 
transfer coefficients at the boundaries. Furthermore, no variation with time or 
temperature is considered in the material properties. 

(b) Numerical Methods 

Most phase change heat transfer problems, particularly those that are 
multi-dimensional, have been solved using numerical techniques. Numerical 
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methods can be divided into two groups, based on the choice of dependent vari- 
ables used: 

1. The temperature is the only dependent variable, and energy conservation 
equations are written separately in the solid region and in the liquid region. 

2. The enthalpy is used as a dependent variable along with the temperature, 
and energy equation is applied over the whole domain covering both solid 
liquid phases. 

Koeger and Ostrach (1974), Springer and Olson (1963), Gooding and Khader 
(1976) , and Lajaridis (1970) used the first approach to obtain numerical solu- 
tions. The second approach, termed the enthalpy method, was used by Sham- 
sunder and Sparrow (1975), Meyer (1970), Dusinberre (1949), and Griggs et 
al. 


Experimental investigation of phase change problems is important in order 
to check the validity of various analytic and numerical models. Thus far, however, 
surprisingly little experimental work has been done and only a limited number 
of experimental studies are available in the literature. Richter (1978) tested a 
thermal energy storage material contained in an annulus bounded by a heat 
pipe. White et al. (1977) performed experiments to observe the melting from 
an electrically heated cylinder cylinder embedded in a phase change material. 
Further experimental work was performed by Thomas et al. (1963), Boger et al. 
(1967), Jiji et al. (1970), and Bailey et al. (1971). 

In the present thesis, first method of numerical method was used for cylin- 
drical annulus geometry to determine the temperature distribution and interface 
location of a PCM using explicit and fully implicit finite-difference scheme. It 
is the purpose of this chapter not only to calculate temperatures and interface 
location but also to evaluate the above mentioned numerical schemes in order to 
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find an efficient finite-different scheme for the phase change problem 

2.2 Fomulation of the Problem 

In the present investigation, attention is given to the following type of 
problem: A given material with a distinct fusion temperature is contained be- 
tween two concentric cylinders. Initially the temperature distribution must be 
known at all points inside the material. Heat is then transfered in or transfered 
out of the material, changing its phase from solid to liquid or vice versa. The 
heat flow in the material is assumed to be such that the fusion front resulting 
from the phase change move towards the outer or inner radial surfaces. Also the 
boundary conditions have to be specified in the form of temperature distribu- 
tions or the heat flow rates along the boundaries of the cylinder and the overall 
heat transfer coefficients at the boundary surfaces. In general these boundary 
conditions may depend on time and on the axial, radial and angular coordi- 
nates. However, in the present treatment only those processes are independent 
of the angular coordinates and consequently where the temperature distribution 
is axially symmetric at all times. Typically in these problems temperature is 
the basic unknown. Once the temperature distribution is known throughout the 
material it is relatively easy to determine other principal variables, such as the 
heat transfer rates and the position of the solid-liquid interface. 

The followings are the assumptions which have been made for the present 
solution method. 1). Phase change meterial(PCM) inside two cylinders is ho- 
mogeneous. 2). Physical properties of PCM are taken to be uniform in the same 
phase. 3). PCM inside two cylinders has a distinct fusion temperature. 4). The 
density of PCM is taken to be same in both solid and liquid phase. 3). The 
natural convection effects are negligible. 


23 


Considering the above assumption, the governing equation for two dimen- 
sional ajrisymmetrical phase change heat-conduction problem is given in the di- 
mensionless form: 


90 _ [9*0 1 90 90 

dr ~ 1 [dR 2 + 1 + R dR + dZ\ 


in Ri<R<E(Z,r ) or E(Z,r)<R<Ro 
and Z j Z Z e 
for r > 0 

with associated boundary conditions and initial conditions: 

1) For r > 0 

|ii = /,(2,x) or 

1®« = /2(Z,t) and B = MZ, t) 


( 2 . 2 ) 


(2.2.o) 


at R = Ri or R = R 0 and Z / < Z < Z e 


and 


90 . „ . 

, az =,,(B ’ T) or 

Qoo = 92 (R,t) and H = gs(R,r) 
at Z — Z f or Z — Z e and R{ R ^ R 0 


2) For r = 0 


Q = fo(R,Z) and 
E = go(R, Z) 


( 2 . 2 . 6 ) 


(2.2. c) 
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in Ri < R < R a and Z f < Z < Z e 

and with the solid-liquid interface condition as 

{ &i = Q, = 0 m = 0 and 

dE _ [00, Jfc,00/1 r (2.2. d) 

6t ~ [6R k t dR J [ + \6Z ) 

at R = E(Z, r) for r > 0 



2.3.1 Finite Differencing of Time Derivatives 

Two different schemes were used to evaluate and compare the results in 
order to estimate the accuracy's and efficiencies (computer storage and CPU time 
requirements) of each scheme. 


(a) Explicit Scheme 

To obtain the stable time dependant two dimensional axisymmetrical heat 
equation, the Fourier stability analysis was performed to yield 


8 MPJZ^dW + dZ^) 

/316dX t +16dZ > + 32<Ut*dZ*+dlPdZ t ' 1 

The derivation of the Equation 2.3 is given in Appendix A. The above equation 
is valid for the axisymmetric two dimensional transient heat conduction problem 
providing periodic boundary conditions are imposed. Since, generally, very small 
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time step size is required to satisfy the Equation 2.3, prohibitively a large number 
of time steps are required for the problems that are to be performed over a large 
period of time. For this reason, an implicit was also investigated to estimate trade 
offs between the time step size (computer CPU time) and degree of accuracy. 

(b) Implicit Scheme 

An implicit finite difference scheme is derived for the problem. Unfortu- 
nately, the resulting system of linear algebraic equations is no longer tridiagonal 
and Thomas algorithm can not be applied because of the five unknowns. The sys- 
tem of equations was solved directly using Gauss Jordan method with maximum 
column pivoting at each time level. 

2.3.2 Finite Differencing of Space Derivatives 

Phase change material of an axisymmetrical geometry is discretized into 
finite difference grids for numerical calculation as in Figure 2.2. Because of the 
discontinuity of the temperature distribution across the solid-liquid interface, 
The computation domain should be subdivided into different regions namely 
inner region, interface region, and boundary region. The choice of the finite 
differencing, therfore, has to be different depending upon grid location. 

(a) Inner Region 

This is the region away from the interface and boundaries, therefore the 
points used in any one of the following difference formulas are not separated by 
the interface or the boundary (Figure 2.3). This region also includes the region 
near the interface and boundary provided that all of the grid points used in 
any one expression are in the same phase. Two or three-point central difference 
formula are used which has second-order accuracy: 
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dQ ^ Qj,j+i ~ Qt,j— i 

dR ~ 2AR 

^ 1 ~ 20 t ,j + Qi,y-i 

5il 2 “ (Ail) 2 

d 2 ® ^ ®*+i,i ~ 2Q»,j + 

0Z 2 “ (AZ) 2 


(2.4) 

(2.5) 

( 2 . 6 ) 


(b) Interface Region 

This region is defined for those three grid points in either radial or axial 
direction: one containing the interface (interface node) and the other two points 
(semi-interface node) directly adjacent to it (Figure 2.4). For these points, pre- 
vious three-point central formulas can not be used because of the discontinuity 
of the temperature at the interface. 

Either three-point forward differencing or backward differencing has been 
used for the grid points containing the interface according to the phase of the 
points relative to the interface. 

The backward difference formula is used when the point lies between the 
interface and the inner radius: 


dR A R 

d 2 e 2 




1 




(2.7) 

( 2 . 8 ) 


OR 2 “ (A R) 2 (2 ' 2 

and the forward differencing is used when the point lies between the outer radius 
and the interface: 


dQ^ JL 

dR A R 

d 2 e 2 


+ 2 ®i,i+i “ 2 €> ‘ ,i+2 


(2.9) 


1 1 

dR 2 ~ (A R) 2 L2 0<li ” ® i,j+1 + 2® i,i+2 


( 2 . 10 ) 
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In Equations 2.7 to 2.10 for the interface grid point used by Springer and 
Olson (1962), temperature of the interface was not included which also means 
the first interface boundary condition of Equation 2.2. d was not utilized. Since 
present author experienced that using Equations 2.7 to 2.10 resulted in numerical 
instabilities in temperatures near the interface, modification has been made to 
remedy that defect. Instead of using three point differencing which excludes 
the interface temperature, two-point (linear) interpolation between the interface 
and the node of same phase adjacent to the interface node is applied as follows 
(Figure 2.5). 


a© 

dR 
d 2 Q 
dR 2 ~ 


AR + 6R 
0 


( 2 . 11 ) 

(2.12) 


Also, for the interface node points near the boundary where the point lies 
between the interface line and the boundary point (Figure 2.6), the same two- 
point (linear) interpolation can be used by replacing A.R-M-R with A Ej. 


6Q 

dR A Ej 


6 2 Q 
8R 2 


a o 


(2.13) 

(2.14) 


where A Ej is the radial distance between the interface and either inner or outer 
boundary. 

For the grid points directly adjacent to the interface points, three-point 
unequal-spacing centered difference formulas are used (Figure 2.7). The differ- 
ence formula is written as 


de ^ 
dR ~ 



- Q%,j - 1 



(2.15) 
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Pe 2 [ 8ij-i _ e,j ' 

dR 2 (AH) 2 [2 + r* 1+r* 


(2.16) 


for the point lies between the interface and the inner radius (Figure 2.7.a). And 
another formula is written as 


de 1 
dR ~ ~ AR 



8*9 s 2 Qj,j-i _ &i,i • 
SR 2 (Al ?) 2 L 2 — r" 1 - r- 


(2.17) 

(2.18) 


for the point lies between the interface and the outer radius (Figure 2.7.b). 

To find the axial derivatives in this region, central difference formula, Equa- 
tion 2.6 is used for point when the two axially neighboring points to the left and 
right of it, are in the same phase. If the points are in different phases, the same 
linear interpolation of Equation 2.33 is used as the radial case. 


(c) Boundary Region 

This region contains the points on the boundaries (Figure 2.8). There- 
fore, the boundary points can be eliminated from the calculation matrix for the 
problems with boundary conditions of the first kind (boundary temperatures are 
prescribed). For the boundary conditions of the second kind (heat fluxes are 
prescribed) and the third kind (ambient temperatures are specified), the tem- 
perature derivatives normal to the boundarys have to be expressed in terms of 
the boundary conditions. 

The first and second derivatives in R direction at the inner and outer radial 
boundaries are given as 


dQ * _ 0 . Un b 

dR ~ Wi ' j r itj k 


(2.19) 


i+1 


the finite 
ier radi 




36 


d 2 ® ~ \ ®t,j±I ®t,j , q, r in kf 

dR 2 “A R [ AR Wi ' j T itj k 

where Q* is the dimentionless heat transfer defined as 


QW 


QjjAQ 

2ir kf A z Tf 


( 2 . 20 ) 


( 2 . 21 ) 


If the interface separates the boundary point from the adjoining one, Equa- 
tion 2.20 can not be applied, in which case the first term on the right-hand side 
of the equation has to be modified. At the inner and outer radial boundaries the 
second derivatives with respect to R are 


d 2 Q 

dR 2 


AEj 


AEj 


4 - Q* . — -L 


( 2 . 22 ) 




Where AEj is the radial distance between the interface and, the inner (j = 1) 
or radial (j = N) boundary 


It should be noted that Equation 2.22 goes to infinity as the interface appoaches 
the inner or outer boundary. Therefore, their computation has to be started 
with the interface already advances enough distance from the boundary to avoid 
this possibility. 

Similar equations are also used to calculate the first end second derivatives 
in the axial direction at the front and the end surfaces of the cylinder: 


d® jjg _q* kf AZ 

dR - r itj k AR 

^ 2 ® 2* £ ®»,j±l ~ ®«,J , Q* r »n kf A Z 

dR 2 ~ AR [ AR * Wi ' j r itj k A 


(2.23) 

(2.24) 


and 
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a 2 €> g 2 

dR 2 A Ej 


®».j , kf A Z 

AEj ViJ rjj Jfc Ai2 


(2.25) 


for the situation when the interface separates the front or end boundary from 
the adjacent point. 

To use the same temperature derivatives written in terms of the heat- 
flow rate Q * for the convection boundary condition, Q* is determined from the 
following relationship: 


Qi,j ~ U i,j AZ r .^ (QooiJ ©t,j) 

where the dimensionless overall heat transfer coefficient U’ is given by 


(2.26) 


it* _ Ujd r * n 

iJ kt 


(2.27) 


(d) Interface Motion 

So far, The derivatives of the temperatures at the node point has been de- 
rived to transform the governing equation (Equation 2.2) into difference form. To 
transform The heat-balance equation (Equation 2.2. d), temperature derivatives 
at the interface have to be determined. At the interface, forward and backward 
interpolation give 


and 



dR A R 




(2.28) 


(2.29) 
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Near the inner or outer radial boundary where only one node-temperature 
is available in any one phase-side (Figure 2.5), the above two equations are to 
be replaced by Equations 2.31 and 2.32. 

At the interior mesh point the slope of the interface in Equation 2.9 is 
determined by the centered three-point interpolation 


9E ^ Ej+i — Ej-i 

dZ ~ 2AZ 


(2.30) 


At those boundaries where i — 1 and i = M the expression for the slope is 


8E ^ ±Ej+i Ej 

dZ A Z 


(2.31) 


2.4 Results and Discussion 

Since there is no analytical solution in two dimensional phase-change prob- 
lem available, the present general method is applied to one dimensional problem 
and the results are compared with approximate analytic solutions. Also a para- 
metric investigation was performed by using different number of radial nodes 
and/or numerical scheme in order to examine grid-dependency of the solution 
and trade-off between cpu time of the computer and solution accuracy. 

Table II shows results of the explicit scheme for the ice formation problem 
of London and Seban. Tabulated values are the elapsed time for solid/liquid 
interface to reach the specified radial positions, namely, J2=1.0, 1.25, 1.5, and 
1.75 for the different number of nodes used (Table Ha). Also required CPU 
time for the solid/liquid interface to reach the outer boundary, r c is tabulated 
in Table lib. Figure 2.9 shows the water/ice interface location in freezing water 
obtained from explicit scheme with 4 radial nodes (solid line). The interface 


39 


location was calculated from the solutions given by London and Seban and was 
also plotted in Figure 2.9 (symbols). Similarly, the water/ice interface locations 
obtained from the explicit scheme with 6, 9, 12, 15 and 18 nodes are plotted in 
figures 2.10 to 2.14 along with analytic solution. 

From Table II and Figures 2.9 to 2.14, the numerical results with more than 
9 radial nodes are shown to be in good agreement with the analytical results of 
London and Seban. The result with 18 radial nodes (Figure 2.14) is not deviated 
far away from the result with 9 radial nodes (Figure 2.11) in Table Ha, while 
required CPU time using 18 nodes are about 8 times larger than the CPU time 
using 9 nodes. 

Table III shows the results of the implicit scheme for the same problem 
discussed above for which the previous explicit scheme has been tested. In this 
table, the same parameters are tabulated as in Table II using the implicit scheme. 
Figure 2.15 shows the water/ice interface location in freezing water obtained from 
implicit scheme with 4 radial nodes (solid line). The interface location was also 
plotted from the solutions given by London and Seban (symbols). Again, the 
water/ice interface locations obtained from the implicit scheme with 6, 9, 12, 15 
and 18 nodes are plotted in figures 2.16 to 2.20 along with analytic solution. 

From Table III and Figures 2.15 to 2.20, it is observed that the implicit 
scheme is somewhat dissipative as the solid/liquid interface advances to outer 
radius for the results with 9 nodes or more. The above diviation can be attributed 
to the growing round off error as the number of the nodes increased due to a 
large number of dependent arithmatic operations. Since two dimentional nature 
of this scheme, the efficiency of the tridiagonal algorithm has been lost and 
the resulting system of algebraic equation had to be solved at each time step, 
one can see, in Table Illb, considerably larger CPU time has been required 
to complete the calculation. Figure 2.16 shows the result with 6 radial nodes 
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is in good agreement with analytic result, however, since the results for other 
nodes show that the implicit scheme is in error, choice of using 6 node implicit 
scheme was not made due to uncertainty for other applications. Now comparing 
the results obtained using the two numerical schemes, explicit and implicit, the 
explicit scheme with 9 radial nodes has been selected to be used in analizing 
the behavior of the phase change material of the thermal energy storage system 
of the receiver. The results of TES analysis will be presented in Chapter IV 
and comparison has been also made with the case without phase change energy 
storage system. 
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Implicit scheme 4 nodes 6 nodes 9 nodes 12 nodes 15 nodes 18 nodes 

Elapsed CPU time, sec 0.5 2.68 19.13 78.46 234.52 577.95 
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Figure 2.9 Comparison of the interface location with the analytic 
result subject to constant temperature: explicit scheme with 4 radial 
nodes. 
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Figure 2.10 Comparison of the interface location with the analytic 
result subject to const amt temperature: explicit scheme with 6 radial 
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Figure 2.1 1 Comparison of the interface location with the analytic 
result subject to constant temperature: explicit scheme with 9 radial 
nodes. 
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Figure 2.12 Comparison of the interface location with the analytic 
result subject to constant temperature: explicit scheme with 12 radial 
nodes. 
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Figure 2.13 Comparison of the interface location with the analytic 
result subject to constant temperature: explicit scheme with 15 radial 







Figure 2.14 Comparison of the interface location with the analytic 
result subject to constant temperature: explicit scheme with 18 radial 
nodes. 
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Figure 2.15 Comparison of the interface location with the analytic 
result subject to constant temperature: implicit scheme with 4 radial 
nodes. 
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Figure 2.16 Comparison of the interface location with the analytic 
result subject to constant temperature: implicit scheme with 6 radial 
nodes. 
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Figure 2.17 Comparison of the interface location with the analytic 
result subject to constant temperature: implicit scheme with 9 radial 
nodes. 
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Figure 2.18 Comparison of the interface location with the ana- 
lytic result subject to constant temperature: implicit scheme with 12 
radial nodes. 
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Figure 2.19 Comparison of the interface location with the ana- 
lytic result subject to constant temperature: implicit scheme with 15 
radial nodes. 
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Figure 2.20 Comparison of the interface location with the ana- 
lytic result subject to constant temperature: implicit scheme with 18 
radial nodes. 
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CHAPTER III 

NATURAL CONVECTION ANALYSIS IN ENCLOSURES 

3.1 Introduction 

Natural convection in open cavity has received growing attention because 
of its important role in solar thermal central receiver system as well as other engi- 
neering systems. Minimizing and predicting energy losses by natural convection 
are major design criteria in such solar systems. During the ground testing of 
the space solar dynamic system, convection heat transfer due to bouyancy force 
occurs inside the cavity receiver between the inner surface and the surrounding 
media i.e., air. Accurate calculation of convection loss through the apart ure is 
need for proper prediction of the receiver heat transfer and further the overall 
performance of the system. 

A review of literture concerned with central receivers has shown that very 
little sophistication is currently being used to analyze convective losses. Com- 
plexity of the problem has necessitated making simplifying assumptions in order 
to solve the problem. However, reviewing the results of these analysis with other 
thermal analyses, it was found that the assumptions axe oversimplified. There- 
fore, more rigorous analytical modeling is necessary for the fluid motion and 
heat transfer and detailed understanding of thermal stratification is needed to 
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be developed. Comparison of the results for receiver convection loss, will be pre- 
sented in Chapter 4. A mathematical model has been developed for predicting 
heat transfer by free convection in a closed container and the results are finally 
compared with available experimental and numerical data. 

This chapter describes an effort to analize natural convection inside enclo- 
sures. Earlier methods applied to convective heating analyses are reviewed. Free 
convection analysis for vertical flat plate has been used in many internal flow 
analysis as a important extreme case. Schmidt and Beckmann (1930) performed 
the first experimental and theoretical work concerning the free-convection flow of 
air subject to the gravitational force about a vertical flat plate, which is consid- 
ered as the most complete treatment of this subject. Eckert (1948) has further 
verified and extended the experimental results of Schmidt and Beckmann, and 
Schuh (1948) has extended the numerical calculations by computing the veloc- 
ity and temperature distributions for several Prandtl numbers. Ostrach (1951) 
obtained the same results as Schmidt and Beckmann via more general approach 
that demonstrates the significance of all the important parameters and assump- 
tions associated with the free-convection flow phenomenon and indicates the 
quantitative limitations of the theory. More information on the free-convection 
flow was obtained from Ostrach’s work for Prandtl numbers corresponding to 
those of liquid metals, gases, liquids, and very viscous fluids. Further analysis 
has been performed by Sparrow (1955) for non-uniform boundary conditions and 
by Siegel for transient free-convection problem utilizing the Karman-Pohlhausen 
method. 

Among some earlier work for flow inside a circular cylinder, Lighthill (1953) 
used analytical techniques in solving the problem of cooling turbine blades. He 
employed an integral- momentum, boundary layer analysis to predict the flow and 
heat transfer due to free convection in heated vertical tubes. Also in his work, 


57 


three different regimes for both laminar and turbulent flow were characterized, 
namely; boundary layer regime, similarity regime and non- similarity regime and 
equation for nusselt number has been derived for each of the flow regimes. For 
turbulent flow, he derive the following simple relation for the boundary layer 
regime by extrapolating previous experimental results. 

Nu r = O.U Ra} r /S (3.1) 

Stratification in cryorgenic propellant storage tanks has stimulated a number of 
papers on the subject of convection in large enclosures where the boundary layer 
flow does not occupy a major portion of the enclosure volume. Clark (1965) 
suggested the subdivision of the flow field into three regions. In his model, a 
flat-plate, natural convection, boundary-layer flow is assumed along the heated 
wall. The flow is discharged into an upper stratified region, and the bulk tem- 
perature below this region is assumed constant. However he noted that one of 
the restrictions of this approach is that the assumption of a constant bulk tem- 
perature is not appropriate except during short start-up time. Later, Evans et 
al. (1968) presented an analytical model using similar divisions fot the flow field 
as suggested by Clark. Hess and Miller (1979) analyzed axisymmetric natural 
convection flows experimently using laser Doppler velodmeter and also devel- 
oped a numerical model with the aid from their experimental data. They also 
noted that both Lighthill and Ostrach did not consider the interaction between 
the boundary layer and the core, which is necessary condition for a complete 
enclosure. 

Schwind and Vliet (1964) studied the natural convection and stratification 
phenomenon of fluids in a vessel with heated sidewalls using the schlieren and 
shadowgraph techniques. They also confirmed the above deficiency noted by 
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Clark in using constant bulk temperature.Their finding was a result of studying 
analogy between pressure gradient in the forced convection and bulk temperature 
gradient in the natural convection. Tatom et al. (1963) performed experimental 
tests under various environmental conditions to study thermal stratification of 
liquid hydrogen in rocket propellant tanks. Hiddink et al. (1976) performed an 
experimental work to develop a better understanding of the mechanism of heat 
transfer within heated liquid in a closed container. This analysis was initiated 
from a special interest in heat sterilization of liquid food in cans or glass jars. 

3.2 Numerical Modeling of Enclosure Flow 

In this analysis, the following assumptions have been made: 

1. The Boussinesq approximation is applicable. 

2. The effect of the wall curvature is negligible. 

3. Boundary layer approximations axe applicable to the flow near the vertical 
walls. 

4. The radial distributions of the vertical velocity and the temperature in 
the boundary layer can be approximated by functions containing four pa- 
rameters (U m , So, T w , and T c ). These parameters are function of vertical 
position and time. 

Reviewing the previous work for the enclosure free convection analysis, 
it has been found that four regions can be identified, namely; 1) the wall, 2) 
boundary layer region, 3) mixing region, and 4) core region for the case of the 
side wall heating. Also another region (low-mixing region) can be added for the 
case of side and bottom heating in addition to the above four regions. In this 
chapter, only side- wall heating is considered. 


3.2.1 Wall of the Enclosure 
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In most practical applications, the ratio of the length of the wall to the 
wall thickness ( L/b ) is very large, and the wall can be treated as a fin. The 
external side of the fin is insulated while the internal side discharges heat to the 
fluid at a rate determined by the local temperature difference between the wall 
and the core of the fluid. As energy balance shown in Figure 3.1 , part of the 
energy from the source is used to heat up the wall and the other part heats up 
the fluid, it is assumed that the characteristic time of heating the wall is much 
smaller than the characteristic time of heating the fluid. 

The governing equation will then be: 


dTre __ a w h (m rr \ 

dt ~ w dz> k w b [ w e) 


(3.2) 


Equation 3.2 is non-dimensionalized using the following parameters: 


t = 



z = LZ 


The final form of the energy equation in the wall is then obtained: 


dT w _ 1 a w d 2 T w 
"dr = A?a 7 dZ 2 


-M r h(T„ -Tc) 


(3.3) 


where, 

L J - r T 

-Mr “ r 7 

tty AJ|u u 


To solve Equation 3.3, the heat transfer coefficient must be known in ad- 
vance. As will be shown in the next section, h can be expressed in terms of the 
boundary layer thickness 6q. In the first time step it is assumed that h = 0. 
The solution to Equation 3.3 is obtained with Thomas algorithm implicit finite 
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difference scheme. The wall equation can then be expressed in finite difference 
form: 


T^ 1 - Tg. 1 a w T2+\ - 2 r£ +1 + 

At i4 2 a/ (Az*) 2 

-X r k"(C, +I -2^) 


(3.4) 


The boundary conditions to Equation 3.4 are: 


dT u 


= 0 at Z = 0 


dZ 

T w = T p at Z = 1 


While initial condition is: 


Tu, = 2/ at t = 0 


The Thomas algorithm is used to invert the tri-diagonal matrix that results from 
expressing Equation 3.4 at every node t. 

3.2.2 The Boundary Layer Region 

The new wall temperature is used to update the buoyant force, and we can 
now proceed to evaluate the boundary layer thickness, So, and the distribution 
of maximum velocity, U m . These governing equations are the integral form of 
the momentum and energy boundary layer equations. The governing equations 
can then be expressed as: 



(3.5) 
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where it has been assumed that 8? = & o • The effect of the dynamic pressure will 
be neglected in the present analysis, although this effect can be very important 
near the top and bottom walls of the enclosure. It is through this dynamic 
pressure that the top and bottom influence the flow field, imposing a zero velocity 
at both ends. 

In order to solve the set of coupled Equations 3.4 and 3.5, we must have 
suitable expressions for the radial variation of the temperature and axial velocity 
inside the boundary layer. These profiles constitute an a priori knowledge of the 
solution, which is obtained from experimental data and/or another solution. 
The unknown parameters then become: the boundary layer thickness (£o)» the 
maximum velocity at any height ( U m ), the wall temperature (T w ), and while 
temperature in the core (T c ). 

The velocity profiles must satisfy the following boundary conditions: 

(1) U = 0 at the wall 

(2) U = 0 at So 

The temperature profiles must satisfy the following boundary conditions: 

(1) T = T W at r = 0 

(2) T = T C at r = 8 0 

In addition, these expression for the velocity and the temperature must 
closely follow the shape of the distributions obtained experimentally. 

The velocity profile can be assumed as following form: 



(i) +C! (i) +Cl (i) +c, (i) 


(3.7) 
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and the temperature profile can be approximated by the following form: 


'Lzh = ( i_JlY 

W -T c V *»/ 


T f - 

T w 


(3.8) 


where a and Ci are usually determined from experiment. 

To generalize the governing equations, we introduce the functions that 
represent the velocity and the temperature profile parameters: 


S' tl 

II 

(3.9) 

T(z,r,t) — T e (z,t) 
2^(2, f) “ T c (2,/) 

(3.10) 


where r 1 = t/Sq. Also following integrals and shear stress 5 are introduced: 

h = f f(r')dr\ I,= [ f* (r')dr', Is = [ g(r') ir’ , 

Jo Jo Jo 

h = £ f(r'Mr')Jr‘, S = 


a/(r') 

dr' 


(3.11) 


I r' = l lr'=0 

the Equations 3.5 and 3.6 are now non-dimensionalized using the following pa- 
rameters: 


t' = ^, U' = 


r T U 


’ * =1’ 


r T a f 

Ty) Tp _$0 A T 

(Ty-T/)’ 71 rx’ 


(3.12) 


Resulting momentum and energy equations axe given as 

2 9 UU , , dg IsVZvth, 
-TT +IiU ^Tt + — ~a? 

+ 2ffi££ = s Pt0: + l u.P,,*/, 


(3.13) 
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and 


V 2 dT c 


, A17 /rr, - dv jjt dO 

+ T (^a? + 


<?z' 


„ du m 

+e ^ m 


f) 


+ 


W ar e 

a at az' 


(3.14) 


= a0 


subject to the following boundary and initial condition: 


U' m = 77 = 0 at z' = 0 

= 77 = 0 at r = 0 

Equation 3.13 and 3.14 are a coupled set of hyperbolic equations for the 
unknowns U ^ and 77 as functions of z' and r. 

In Equation 3.14, the value of 9 must be provided by the solution of the 
wall equation. This set of coupled, non-linear, partial differential equations are 
solved using Newton-Raphson iteration method. 

3.2.3 The Mixing Region 

The mixing region is a control volume in the uppermost part of the con- 
tainer, where the boundary layer discharges its energy. It is bounded by the 
boundary layers and has a depth E. The mixing region is considered to be 
important region through which the flow coming from the boundary region is 
entering into the core region, since unlike flat plate case, the velocity and the 
temperature of the boundary layer flow will affect those of the flow outside the 
boundary layer (core region) via mixing region in the enclosure problem like 
this. Therefore accurate prediction of the flow characteristics in that region 
would improve the solution for the whole heating process. 

An energy balance in this control volume (Figure 3.1) provides the average 
mixing temperature T m . This temperature serves as the boundary condition to 
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the core. It is assumed to be in the middle of the mixing region. The energy 
balance states that: the energy convected by the boundary layer at H plus the 
heat convected into the fluid between the top and H equals the energy stored in 
the mixing region plus the energy that leaves the mixing region at a velocity U c - 
That is: 


o rL 

' p/UCpfTf 27r(rj — r)dr + 2tt r? / q v dz 

o g-L-H J L-H 

= p, Op, T(r T - S 0 + p, C, 2*^ f* (T u - T,)( r T - r) ir (3.15) 

r^T 

- / pfU c CpfTm2n(r T -r)dr 

J6 n 


z=L-H 


where the heat flux q w is obtained by 


g. = ^(r„ - r e ) 


The initial condition to Equation 3.9 is: 


Tm = Tj at t = 0 

The solution of Equation 3.9 provides Tm, which is assumed to be the tempera- 
ture in the center of the mixing region (z — L — H/ 2). It is still not clear what 
is the best method of estimating the temperature profile above the middle of the 
mixing region. It is reported that parabolic temperature profile for this region 
under-estmates the temperatures near the top which results in a larger driving 
force than actual situation. Consequently, instead of parabolic profile, following 
power function has been used for the temperature profile with prescribed T m at 
the mixing region. 
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y = a/(x 2 + b) 

A special form of tHis curve is called ‘Witch of Agnesi’ when a = 6 3 / 2 , and 
these curves automatically satisfy the boundary condition at the top (insulated). 
In Figure 3.3, excellant agreement has been obtained with this profile for the 
same test case of Hess (1982). 

3.2.4 The Core Region 

The temperature of the core is assumed to be constant in the radial direc- 
tion, which is close to the behavior shown by the numerical solution. It is allowed 
to vary with the axial coordinate and with time. The equation that defines the 
temperature in the core is: 


dT t 


6T e 8 2 T e 


dt +Ue dz a/ dz* 


(3.16) 


where the velocity of the core can be closely expressed by: 


U c = 


-2 


( r T - 6 0 ) 2 


o 

/ U(r T -r ) 
Jo 


dr 


z=L—H 


Neglecting the effect of curvature, the core velocity can be rewritten as: 


2 f 6 ° 

U c = — jr* / Udr 
■*w> Jo 


z—L—H 


and using the definitions in Equation 3.11 it can be expressed as: 


U c = -2 h U m 


Vm&x 


max 


(3.17) 
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Equations 3.17 and 3.18 are non-dimensionalized using Equation 3.12. The re 
suiting expressions are: 


6T C U' c dT c _ 1 8 2 T C 
dr + A dz' A 2 dz' 2 


(3.18) 


and 


K = -2 hU' m 


Wx 


The boundary conditions to the energy equation are: 


(3.19) 


T c — Tm at z — 1 

The solution to Equation 3.19 is obtained with an implicit finite difference 
scheme. The method proceeds by updating the wall temperature which, in com- 
bination with the core temperature of the previous time, provides the driving 
force, 6, necessary to solve the parameters of the boundary layer. 


3.3 Results and discussion 

The results obtained with the integral analysis with revised temperature 
distribution in the mixing region are varified by comparing with those from Hess 
(1982) and from Miller (1977) who solves the same problem using full Navier- 
Stokes equation. Comparison of the velocity profiles with the result of Hess for 
the sample cases are presented in Figures 3.2 to 3.5 for two different Rayleigh 
numbers. In the Figures, solid lines with symbols represent present results , and 
dashed lines and symbols represent the results from Miller and Hess respectively. 
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Figure 3.2a shows the variation of the velocities vs. enclosure height at 
three different times, namely, 60, 180 and 360 seconds for Ra=3.73 x 10 8 . Also, 
Figure 3.2b shows similar results for Ra=7.46 x 10®. From these figures, it is 
noted that the present results compare fairly well with the numerical solution 
and matchs to the Miller’s numerical results better than Hess’s integral results 
particularly for later time (<=360 sec.). 

Figure 3.3a shows the variation of the core temperatures vs. enclosure 
height at the three different times mentioned above for Ra=3.73 x 10®. Also, 
Figure 3.3b shows similar results for Ra=7.46 x 10®. From these figures, it is 
noted that, with the use of power function for temperature profile in the mixing 
region, excellent agreement has been obtained with the Miller’s results at all 
three times. 

Again, Figure 3.4b shows the variation of the wall temperature vs. en- 
closure height at the three different times mentioned above for Ra=3.73 x 10®. 
While Figure 3.4b shows similar results for Ra=7.46 x 10®. The present integral 
results compare well with the Miller’s solution except in the early time. That 
can be attributed to the fact pointed by earlier investigators that the boundary 
layer assumption is not valid at early times. Accordingly, this analysis can be 
appropriate for the solar receiver analysis which involves cylic behavior which 
requires substantial period of time to achieve steady results. 

In Figure 3.5a , the variation of the local Nusselt number vs. enclosure 
height at four different times (60, 120, 300, 420 sec.) are plotted for Ra=3.73 x 
10®. Also, similar results is shown in Figure 3.5b for Ra=7.46 x 10®. The present 
solution is compared with the numerical results and very good agreement was 
obtained. The following can be noticed from this figure: 

1. ) There is very little variation of local Nusselt number with time. 

2. ) The variation is almost linear which implies that h is almost a constant. 
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From the previous figures, it verifies the current analysis compares fairly 
well with the earlier numerical results and suitable for the receiver analysis which 
deals with cyclic variation of the parameters. Application of this analysis to the 
free convection analysis inside a receiver cavity will be presented in Chapter IV 
along with results obtained. 
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Figure 3.2 Variation of the velocities vs. enclosure height at three 
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CHAPTER IV 

APPLICATION TO SOLAR DYNAMIC RECEIVER ANALYSIS 

4.1 Introduction 

Numerical technique presented in Chapter II and III has been applied to the 
receiver analyser computer program, HEAP, which was discussed in Chapter I. In 
this chapter, more detailed description of the receiver heat transfer calculations 
of the HEAP code is made followed by the sample simulations of a solar receiver 
under the different operating conditions, namely, one on the orbit and the other 
on the ground. 

4.2 Finite- Difference Receiver Program Methodology 
4.2.1 General Layout 

The receiver under consideration is discretized in space by nodes as shown 
in Figure 4.1. In this figure, each component of the receiver is represented by 
the finite element nodes as follows: 

o Nodes 1 to 10: Working gas of the Brayton cycle, 
o Nodes 11 to 20: Working gas tubing, 

o Nodes 21 to 30: Phase change material of the thermal energy storage 
device, 
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o Nodes 31 to 40: Thermal energy storage containment, 
o Nodes 41 and 42: Aperture plate, 
o Nodes 43 to 45: Back plate, 
o Node 46: Working gas inlet, 
o Node 47: Aperture opening. 

Heat exchange between any two arbitrary nodes (i) and (j) is illustrated 
in Figure 4.2. The heat transfer rate equations are written with all modes of 
heat transfer included. The net energy stored in each node is calculated using 
the first law of thermodynamics, including the energy exchange to and from each 
node as a result of neighboring nodes. 

For the transient solution, explicit numerical technique is used. By this 
technique, the new temperature at time (r + At) is determinded using the old 
temperature vector at time (r) and all other parameters evaluated at time (r). 
to ensure stable computations towards convergence a special formula is used for 
the critical time step (At) based on the methodology stated in reference by 
Razelos (1973). 

4.2.2 Basic Assumptions 

The following assumptions and idealizations were made in the heat transfer 
calculations: 

(1) The physical system under study is treated as a collection of a number of 
isothermal nodes with uniform optical and physical properties throughout 
each node. Nodes are assumed to be of axisymmetric geometry accordingly 
one can make use of the special radiation view factor subroutine. The axis 
of symmetry must coincide with the receiver axis. 

(2) The nodal specific heat, thermal conductivity, and density are assumed 
independent of operating temperature. The properties may be taken at 
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Figure 4.1 Receiver Discretization into Nodes. 
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Figure 4.2 Heat exchange between node (i) and (j). 
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an estimated average working temperature as an approximation. How- 
ever, upon consideration of the natural convection in the receiver cavity, 
density of the air is taken function of the temperature according to the 
ideal gas equation. Also, nodal infrared and solar emittances are assumed 
independent of both directivity and operating temperature. 

(3) The radiation emitted by a node surface is assumed to be diffuse. In 
general, the radiation reflected from a node surface can be diffuse, specular 
or a combination of diffuse and specular parts. In the present model, only 
diffuse reflection was assumed, with a zero specular component. 

(4) The optical properties of a node surface are assumed to be divided into 
two bands in the eletromagnetic spectrum. This “semi gray” assumption 
allows separate calculations of radiation energy exchange in the short wave 
(solar) band and the long wave (infrared) band. Consequently, the radia- 
tion energy in the program is allowed to be absorbed, reflected or emitted 
in the long wave region, and only absorbed or reflected in the short wave 
region. 

(5) The nodal internal energy generation or dissipation by eletrical, chemical or 
nuclear effects is assumed uniform over the node surface and independent 
of the temperature. 

Based on the above assumptions and idealizations, the heat transfer rate equa- 
tions are written and analyzed as described in the following sections. 

4.2.3 Nodal Energy Balance 

For a particular receiver node (i), the energy is exchanged with the neigh- 
boring nodes by many modes of heat transfer such as: 

1) Conduction, Qeond'i 

2) Thermal convection, Q conv ; 
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3) Fluid convection, Qf, 

4) Long wave (infrared) radiation, Qir] 

5) Short wave (solar) radiation, Q t ] 

6) Internal heat generation (or dissipation) by electrical, chemical or nuclear 
effects, Q(i) p - 

The first law of thermodynamics applied to the ‘i’th node gives: 

DEN(i) • V(i) • CP{i) = BES{i) = 

Q(i)lR + Q(i)cond + Q{ 0 eonv + Q(i)f + Q(i)$ + Q(i)p 

where BES(i), DEN(i), F(i) are the net. energy stored, density and volume of 
the ‘i’th node, respectively, and At is the elapsed time increment. 

In the following sections, detailed description will be given for the different 
terms in Equation 4.1. 


(4.1) 


4.2.4 Conduction Hea.t Transfer 

In this section, the conduction heat transfer analysis will be discussed. 
For two neighboring nodes (i) and (jf), exchange heat by conduction across the 
interface area A(i,j), the heat transfer Q(i) C ond is given by 


Q(i)c^ = £ C(i, • [ T(j) - r(0] 


(4.2) 


where 


fit 2 I •A(iiJ) 

- 1 -jsn~sz 


W) + Kli) 


(4.3) 


K is the thermal conductivity, and s is the ‘conduction distance between the 
center of each node and the surface separating the two nodes, t and j' (see 
Figure 4.3a). 


C-2. 
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Since conduction heat transfer inside the thermal energy storage (TES) 
material has been identified to be involved in complex phenomena as discussed 
in Chapter II, the Equation 4.2 simply can not be used between two TES nodes. 
Hereby, conduction heat transfer inside the TES material is treated seperately in 
the TES subroutine and the energy balance described in the previous section, is 
not performed for the TES nodes instead temperature has been set to the fusion 
temperature of the phase change material. However, conduction heat transfer in 
the radial direction need to be calculated between phase change material node 
and either thermal energy storage containment node or gas tube node in order to 
determine the temperatures of the latter nodes (i.e., TES containment nodes or 
gas tube nodes) using Equation 4.1 (see Figure 4.4). It should be noted that the 
distance s used in Equation 4.3 is no longer measured from the center of the node 
for the TES nodes. Thus, for TES nodes, s is defined as the ‘conduction distance 
between the solid-liquid interface (or fusion front) and the surface seperating the 
two nodes, i and j' (see Figure 4.3b). 

The value of C(i,j) could be a fuction of time and/or temperature or 
any other variable. However, the conductance C(i,j ) is assumed constant and 
uniform throughout the node. 

4.2.5 Convection Heat Transfer 

The convection heat transfer across the boundary layer of solid-fluid in- 
terfaces (working gas tubing), can be put in the following form between two 
neighboring nodes (t) and (j) of different phase: 

Qli) cn, = £ ■ [T(j) - T(t)] (4.4) 

j 

where the conductance C(i,j) conv is given by: 
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Figure 4.4 Connection of the heat transfer calculation between 
TES material and Receiver, (a) Energy balance for receiver nodes 
adjacent to thermal energy storage node, (b) Finite difference nodes 
for TES calculation. 
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C(i,j) con v = h(i,j) • A(i,j) (4.5) 

and where h(i,j) is the convective heat transfer coefficient. 

C[i,j)conv is treated in a manner that allows its value to be changed with 
time. This is made simply because the convection coefficient h(i,j) is a dominant 
function of the fluid mass flow rate which may not be constant during the receiver 
operation. 

4.2.6 Fluid Convection 

Consider a node (i) surrounded by two adjacent nodes (s) and (k) with 
a fluid flowing at a rate rhf from (s) to (t) to (fc) as shown in Figure 4.5. An 
energy balance on node (i) taking into account the energy entering at T(s) and 
leaving at T(i) would yield 

0(0/ = *M0 • cm*' ) • [ n> ) - r(i) ] (4.6) 

which will be written as 

<}(<)/ = c(i,.) / .[r(.)- no] (4-7) 


or 


0(0/ = E 0(4. i)/ • i T <j) - r (4)] (4- 8 ) 

3 

where T(j) is the temperature of all neighboring nodes in direct contact with 
node (i). For a single fluid stream crossing node (i), there will be only two sur- 
rounding nodes, one upstream and the other downstream. The flow conductance 



Figure 4.5 Fluid Flow Across Node (i). 
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C(i,j) will then be taken els [m/(z) Opf(i)] if the flow goes from ( j ) to (z), and 
zero if the frow goes from (z‘) to (j). The matrix C(i,j)f will then be asymmet- 
ric such that C(i,j) ^ C(j,i) and all elements are needed els input data. For 
boundELry nodes that present the inlet fluid zone or outlet fluid zone els shown 
in Figure 4.1 the next neighboring node in the closed loop should always be 
considered. 

4.2.7 Infrared. (IR) and Solar Radiation 

The analysids of heat transfer by radiation is divided into two wavelength 
regions: a) an infrared (or long wavelength) region in which infrared proper- 
ties apply; and b) solar (or short wavelength) region in which solar radiation 
properties apply. 

The net long wave (infrared) radiation received by the ‘z’th node from all 
neighboring nodes is given by 

Gw™ = E py) - ry)] (4.9) 

j 

where T is the absolute temperature and C(i,j)jR is the IR radiation conduc- 
tance between nodes (z) and (j). From IR radiation analysis, C(i,j) can be 
written as 


C(i,j)lR = A(i) • <r • £(>) • F(,,j) ■ t (i) 

where a is the Stefan-Boltzman radiation constant, e is the gray emissivity of the 
‘z’th node, and F(i,j) is sometimes called total view factor since it does include 
not only the the nodes z' and j but also all the neiboughboring nodes, effect of 
the neiboring nodes as well. 
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Reflectivity and emmissivity have been refered to as p* and e* respectively 
in the solar band as p and e respectively in the IR band. Thus the net solar 
energy exchange on node (z) from all neighboring nodes is given by 


Q(i)s = £ e*(z) • F*(z, j) • p(jy . Flux(j) + e*(z) • Flux(i) • A(i) (4.10) 


where Flux(i) is the incident solar flux on node (i). 


Detailed analysis for the radiation conductance calculation can be found 
in Lansing (1979). 

4.2.8 Transient Temperature Calculation 

After the elapse of a small time interval At, Equation 4.1 gives the slope 
of temperature array T for each node as an approximation of the temperature 
derivative using finite difference. The explicit scheme is used and the new nodal 
temperature T(i) at (r + At) is calculated using the old temperature t(i) and 
all other properties at the old time (t), i.e., 


(4.1D 

where M(z‘), Cp(i) and KES(i) are the mass, specific heat and energy residue at 
the zth node. 

4.2.9 Matrix Formats 

To facilitate the manipulation process, Equations 4.2, 4.4, 4.8, and 4.9 are 
put in matrix or vector forms. Generally, 


Q(i) = £ C(i,j) [!“(;•) -no] 

i 


(4.12) 
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where n is an exponent that can be either 1 or 4, and C(i,j ) is a general form 
for conductance matrix. Equation 4.12 can be written in the expanded form 

<?(;) = <?(«, i) [r*(i)-r*(i)]+c(<.2) [r*(2)-r*(o] +••• 

or 


Q(i) = c(i, i) r n (i) + c(i, 2) r n (2) + . • • - [c{i, i) + c(i, 2) + c (»\ 3) + • • •] r n (i) 


In the matrix form, Equation 4.12 may also be expressed as 


Q{i)= (C(t,l) C(i, 2) ... -[C(t\l) + C(i,2) + H ...) 


' r n (i) ' 
r n (2) 

T n (i) 


V • / 

(4.13) 


or 


4 ( 0 -£<?(•'. 

i 

where [C] is a modified conductance matrix developed from the original [C] 
matrix with diagonal elements equal to the negative sum of the corresponding 
row. Therefore, for all of the receiver nodes (for i — 1 , 2, 3, . . .), Q{%) can be 
written as the following matrix form: 


[<?] = Pi ■ [t*] 


where [C] is equal to 
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'-[C(l,2)+C(l,3)+-] 

C(2,l) 

C(3,l) 


C(l,2) 

-[C(2,l)+C(2,3)+-] 

C(3,2) 


C(i, 2) 


C(l,3) 

C(2,3) 

-[C(3,l)+C(3,2)+-] 
C(i, 3) 


C(l,i) 

C(2,») 

C(3,») 

-[C(t,l)+C(t,2)+-] 


Based on the above discussion, Equations 4.2, 4.4, 4.8, and 4.9 can be 
written in terms of the modified conductance matrices as 


Q(i)lR = '£C{i,3)lR-T i U) ' 

i 

Q{})cond = ^ ^ j)cond * 2"*(j f) 

_L, _ ► 

Qifyconv = ^ ^ j)conv ’ ) 

j 

3 > 

and the nodal energy balance, Equation 4.1 is written as: 


(4.14) 


*ES(») = Y. ■ T i U) 


+£ 


^(^j)cond + conn + <5(i,i)/ 


r(j) 


(4.15) 


+ Q(0* + Q{*)r 


4.3 Receiver Orbital Simulation with Thermal Energy Storage 


4.3.1 Incorporation of TES Analysis into the HEAP Program 


Phase change thermal energy storage analysis of the Chapter II has been 
incorporated into the HEAP program as a subroutine. In Figure 4.4. a, energy 
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balance for a particular receiver node is depicted in accordance with Equation 4.1 
and for the calculation of the receiver node temperature in main HEAP program, 
the radial TES nodes from j = 1 to N in each column of Figure 4.4b is repre- 
sented by one node whose temperature is set to the fusion temperature. 

Connection of the heat transfer calculation between TES material and 
Receiver is made by the following algorithm for node (i) (see Figure 4.4): 

o Step 1.) With the initial or previous interface location, thermal conduc- 
tance between the nodes (t) and (j) is determined from Equation 4.2 and 
further Qt,eond is determined, (in HEAP program) 
o Step 2.) With the Qt,eond and othe modes of heat transfer to node i, energy 
residue is calculated from the Equation 4.1 and further new temperature 
for the receiver node (i) is obtained, (in HEAP program) 
o Step 3.) After evaluating new temperatures of all the receiver nodes, phase 
change problem is solved for TES nodes (Figure 4.4.b) using new receiver 
temperature as boundary condition, (in TES subroutine) 
o Step 4.) After evaluating new locations of the interface for all nodes, repeat 
step 1 to 3 for next time step. 

o Step 5.) At the completion of the whole cycle (sun + shade time) steps 1 
to 4 are repeated until steady cyclic operation take place. (About 4 cycles 
are needed to reach this conditions) 

4.3.2 Results of Orbital Simulation 

Execution has been performed with the practical environmental condi- 
tions which simulates the actual orbital condition. Using a computer program 
called CAV, incident solar flux distribution in an axi sy mm etrical cavity from a 
paraboloid concentrator is determined for a given solar intensity. In order to 
estimate the effects of the thermal energy storage on the temperature of the 
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various receiver parts, the results with TES and without TES axe compared in 
Figures 4.6 to 4.11. 

Figure 4.6 shows the working fluid temperature distribution vs. time at 
different axial nodes representing inlet, middle and outlet section of the working 
gas with phase change material (PCM) included. In order to estimate the effect 
of PCM on the working fluid temperature, similar results to the above is shown 
in Figure 4.7 with no PCM. Comparing Figures 4.6 and 4.7, the following can 
be noticed: 

1. ) The inlet temperatures are almost same 

2. ) The minimum temperature of the working gas is about 930 °F without 

PCM and is increased about 10 °F in the presence of PCM. 

3. ) The maximum temperature of the working gas is about 1300 °F without 

PCM and is reduced to 1260 °F (40 °F decrease) in the presence of PCM. 

4. ) The temperature fluctuation during a cycle can be as high as 170 °F with- 

out PCM and is reduced to 60 °F due to PCM. 

Also, Figure 4.8 shows the temperature distribution of the working gas heat 
exchanger tubing vs. time at different axial nodes representing the inlet, middle 
and outlet section of the working gas tube with PCM included. The effect of 
PCM on the gas tubing temperature can be observed by comparing this figure 
with Figure 4.9 which shows similar results with no PCM included. Comparing 
Figures 4.8 and 4.9, the following can be noticed: 

1. ) The maximum temperature of the gas tubing is about 1380 °F without 

PCM and is reduced to 1350 °F (30 °F decrease) in the presence of PCM. 

2. ) The minimum temperature of the gas tubing is about 1050 °F without 

PCM and is increased to 1150 °F (100 °F increase) due to PCM. 

3. ) The temperature range can be as high as 230 °F without PCM and is 

reduced to 95 °F with PCM included. 
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4.) The maximum temperature range of the gas tubing occurs at gas inlet 
section as opposed to gas temperature for which maximum temperature 
occurs at gas outlet section. 

Figure 4.10 shows the temperature distribution of the thermal energy stor- 
age (TES) containment vs. time at different axial nodes located at the inlet, 
middle and outlet section of the working gas tube with PCM included. Again, 
the effect of PCM on the TES containment temperature can be observed by 
comparing with Figure 4.11 which shows similar results with no PCM included. 
Comparing Figures 4.10 and 4.11, the following can be noticed: 

1. ) The maximum temperature of the surface temperature is about 1480 °F 

without PCM and is increased to 1750 °F (270 °F increase) in the presence 
of PCM. 

2. ) The minimum temperature of the gas tubing is about 1080 °F without 

PCM and is increased to 1409 °F (329 °F increase) due to PCM. 

3. ) The temperature ranges in the two results are almost same. 

4. ) The minimum temperature of TES containment has been droped to its lim- 

iting temperature which is set purposely equal to fusion temperature of the 
PCM in order to prevent another solidification neax the TES containment. 
Solid-liquid interface location vs. time is plotted in Figure 4.12 at four 
different axial locations, and percent of PCM melted at the middle section of the 
gas tube is shown in Figure 4.13. From the Figures 4.12 and 4.13, the following 
can be noticed: 

1. ) Steady cyclic variation was obtained after first cycle. 

2. ) At the end of shade time most PCM solidifies, however, only about 60 % 

melts at the end of sun time. 

If change of gas mass flow rate is allowed in overall power system design 
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requirements, decreasing the working gas mass flow rate during the sun time is 
suggested in order to increase percent of PCM melted. 

Also, two dimensional motion of the interface vs. axial location is shown 
in Figures 4.14 and 4.15 during sun period and shade period repectively. From 
these figures, the following can be identified: 

1. ) The axial distribution of solid/liquid interface resembles that of incident 

solar flux during the sun time. 

2. ) The PCM solidifies more rapidly during the shade time than it melts during 

the sun time. 

3. ) At the end of shade time, most PCM solidifies completely except the region 

located at the outlet section of gas tube. 

Figure 4.16 shows energy transfered in the receiver vs. time over the cycles. 
The square symbols represents the accumulated solar energy in and circular and 
triangular symbols represents the accumulated and dispersed energy extracted 
by the working gas respectively. The following can be noticed in this figure: 

1. ) During the sun time constant solar energy is coming in and no energy is 

coming in during the shade time. 

2. ) Almost same amount of energy is extracted during sun and shade time. 

In Figure 4.17, variation of receiver efficiency vs. time over the cycles is 
plotted. The receiver efficiency is defined as follows. 


VReeeiver — 


Accumulated Energy Extracted by Working Gas 
Accumulated Solar Energy In 


(4.16) 


From this figure receiver efficiency is approaching about 90 %, and it is increasing 
during the shade time since no energy is coming in and energy is still extracted 
by gas from thermal energy storage. 


93 


Table IV shows a summary of temperature variations of the receiver compo- 
nents shown in Figures 4.6 to 4.11 in order to elucidate the comparison between 
the TES receiver results and the one without TES. The results without TES was 
obtained by replacing the thermal properties of the PCM nodes by the prop- 
erties of TES containment to maintain the same area of heat transfer toward 
the working gas tubing. Therefore one can consider the case without TES as 
the case with sensible heat storage where storage material has the same thermal 
properties as TES containment. 

Since the design objective for the thermal energy storage is to minimize 
time variation in working fluid outlet temperature by providing continuous heat 
into the working gas during the shade time, gas outlet temperature was reviewed 
first. It has been obserbed that using phase change material as a thermal energy 
storage reduces cyclic variation of the working gas outlet temperature down to 
60 °F which is about 65 % reduction from the result obtained in the case without 
TES (170 °F). Consequently, the adventage of latent heat storage system over 
sensible heat storage system discussed in Chapter I can be verified by obtaining 
less temperature variation of the working fluid. However it is observed that 
almost same time averaged temperature was obtained in both cases. 

Similarly, the maximum temperature variation of the TES working gas heat 
exchanger tubing has been reduced about 95 °F which is 50 % reduction from 
the result obtained in the case without TES (230 °F). Maximum temperatures of 
the gas and gas tubing is higher when no TES are used. On the contrary, for the 
TES containment, one can observe reversed behavior. The overall temperature 
of the cavity surface temperature for the phase change TES case is higher than 
that without TES case. This is due to the large heat capacity and higher thermal 
resistance of the phase change material. While due to the relatively small heat 
capacity and thermal resistance, heat transfer rate is higher in the case without 
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TES, resulting lower cavity surface temperature. 

4.4 Receiver Terrestrial Simulation with Natural Convection 

4.4.1 Results of Preliminary Natural Convection Analysis 

As mensioned earlier in Chapter III, preliminary NC analysis has been 
conducted using constant convective heat transfer coefficient from the following 
formula (Equation 3.1): 


Nu r = — = 0.11 Raj/ 3 (4.17) 

k f 

One of the results for closed Brayton cycle (CBC) solar receiver simulations is 
presented in this section. Plotted values in Figures 4.18 and 4.19 axe the results 
of steady-state calculation of the sun period. 

Figure 4.18 shows temperature distribution of the cavity surface vs. axial 
nodes. The solid line represents the results with natural convection and dashed 
line represents the results without natural convection included. These results are 
obtained by steady-state calulation during the sun period, therefore no transient 
behavior for shade period has been shown. Due to natural convection lowest 
TES containment temperature has been dropped about 14 % and highest TES 
containment temperature has been dropped about 7 %. 

Similarly Figure 4.19 shows temperature distribution of the working fluid 
vs. the axial nodes. The working gas temperature has been dropped about 9 % 
due to natural convection heat loss to the ambient air and as a result, receiver 
efficiency has been dropped about 22 %. 

Input values for the evaluation of the h value from Equation 4.16 and 
comparison of the two outputs are tabulated in Table V. It shows the results 
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of the preliminary study of the free convection effects on the receiver using the 
convective heat transfer coeffcient value from the Lighthill’s emperical formula, 
Equation 3.1. The values of parameters in Equation 3.1 which is listed in Ta- 
ble Va, are evaluated at the mean temperature of the ambient air and average 
cavity wall temperature. Resulting output of the receiver simulation are shown 
in Table Vb, and calculations were made only for steady-state during the sun 
time. From Table Vb, due to natural convection, gas outlet temperature drops 
81 Deg. C, which is about 16 kW of heat is lost through the air constantly dur- 
ing the sun time. Consequently overall cavity surface temperature and receiver 
efficiency have been droped also as shown in Table Vb. 

In the following section detailed description on the application of the en- 
closure natural convection problem will be presented along with the results. 
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4.4.2 Incorporation of Natural Convection Analysis into the HEAP Program 

In the previous preliminary study, it is identified that certainly natural 
convection affects the receiver thermal behavior, however, convective heat trans- 
fer coefficient used in this analysis is from the analysis of vertical tube without 
closed upper plate, whereas receiver cavity has the closed upper plate. Therefore 
different flow analysis is necessary for the receiver cavity geometry in order to 
obtain more practical results. 

It was found that air outside the receiver aperture has little influence on 
the air flow inside the cavity. Based on that assumption, the receiver cavity 
can be assumed to be a complete cylindrical enclosure and then the analysis 
discussed in Chapter III can be applied to the receiver free convection analysis. 
Utilizing this analysis, following can be achieved: 

1. ) h will vary with axial direction; 

2. ) h will be time-dependent throughout the cycle; 

3. ) h will be calculated from the analysis that properly characterizes the fluid 

flow inside an enclosure. 

Internal natural convection analysis of the Chapter III has been incorpo- 
rated into the HEAP program as a subroutine. In Figure 4.20, vertically placed 
receiver is shown and cavity is divided into three regions: boundary layer, mix- 
ing and core region as was done in Chapter III. However, some modification has 
been made to the analysis for receiver application as follows: 

(1) The wall region haw been removed from the analysis since heat transfer 
analysis of the cavity wall is taken care by the HEAP program; 

(2) The mixing region has been confined to the upper 10 % of the receiver in 
accordance with the analysis of Clark (1964) and Evans (1968) for constant 
heat flux boundary condition; 
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(3) Exponent of the temperature profile in Equation 3.8 has been changed to 2 
in accordance with the temperature profile of Karman-Pohlhausen method 
used by Siegel (1958) in his vertical flat plate free convection analysis. 
Connection of the heat transfer calculation between the air inside receiver 
cavity and the cavity wall is made by the following algorithm for node ( i ) (Fig- 
ure 4.20): 

o Step 1.) With the initial or previous local convective heat transfer coef- 
ficient, thermal conductance between the node (i) and ( j ) is determined 
from Equation 4.3 and further Q a ,eonv is determined, (in HEAP program) 
o Step 2.) With the Q a ,eonv and othe modes of heat transfer to node t, energy 
residue is calculated from the Equation 4.10 and further new temperature 
for the Receiver node (i) is obtained, (in HEAP program) 
o Step 3.) After evaluating new temperatures of all the receiver nodes, nat- 
ural convection (NC) analysis is performed for the air nodes using new 
receiver temperature as boundary condition, (in NC subroutine) 
o Step 4.) After evaluating new local convective heat transfer coefficients, 
repeat step 1 to 3 for next time step. 

o Step 5.) At the completion of the whole cycle (sun + shade time) steps 1 
to 4 are repeated until steady cyclic operation take place. (About 6 cycles 
are needed to reach this conditions) 

4.3.2 Results of Terrestrial Simulation 

The temperature variations of the air inside receiver core and mixing region 
during the cycles axe plotted in Figures 4.21 to 4.26 with the wall temperature. 
Figures 4.21 to 4.25 shows the core temperature distribution of the air vs. time 
at various axial locations, namely, z/L=.l, .4, .7, .9 and .95. The calaulation 
starts with air temperature inside the cavity equals to the ambient air outside 
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Figure 4.20 Connection of the heat transfer calculation between 
the air and cavity wall, (a) Energy balance for receiver nodes adja- 
cent to the air node, (b) Finite difference nodes for free convection 
calculation. 
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the receiver aperture. Therefore, this is cold staxt-up simulation of the receiver 
for the air. From these figures, the following are shown to be noticed: 

1. ) The core temperature increases rapidly during first and second cycles and 

shows rather steady variation after fifth cycle. 

2. ) It is shown that during the shade time the temperature of the core is higher 

than the average side wall temperature which indicates that the core of the 

air is not as sensitive to the solar heat flux input as the side wall. 

3. ) The mixing temperature is almost identical to the average wall tempera- 

ture. 

Figures 4.21 to 4.25 plotted in one graph in Figure 4.26 in order to view 
the temperature distributions over the entire axial length. From this figure, it 
is noticed that, as we approach to the bottom of the cavity where aperture is 
located, maximum temperature over one cycle is dropped and occurs later time, 
i.e., There is phase lag for the maximum core temperature and is increasing as 
z decreases, however, the amplitude of the maximum temperature is decreasing. 

Also axial variations of the air temperature are shown in Figures 4.27 
to 4.29. Figure 4.27 shows the axial variation of the air temperatures during sun 
time, the solid line represents the temperature distribution at the biginning of 
the sun time and the following dashed lines are the temperature distributions 
at 7 minutes after the previous plot time. From this figure, the follow can be 
noticed: 

1. ) During the beginning period of sun time, temperatures of upper-half por- 

tion of the core increase rapidly but temperatures are still decreasing in 

lower-core region. 

2. ) Temperature range between maximum and minimum is much higher in the 

upper-core than in the lower-core region during the sun time. 
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3.) Maximum temperature at each time step occurs at the top of the core, 

and minimum temperature occurs at the bottom of the core except at the 

beginning of the sun time. 

Similarly, Figure 4.28 shows the axial variation of the air temperatures 
during shade time, the solid line represents the temperature distribution at the 
biginning of the shade time and the following dashed lines axe the temperature 
distributions at 7 minutes after the previous plot time. From this figure, the 
follow can be noticed: 

the following can be identified from the Figure 4.28: 

1. ) During the beginning period of shade time, temperature of the upper-core 

is decreasing, but temperature of the lower-core is still increasing. 

2. ) Temperature range between maximum and minimum is much higher in the 

upper-core than in the lower-core region during the shade time. 

3. ) As opposed to the previous plot during the sun time, maximum tempera- 

ture at each time step occurs at the middle or the bottom of the core, and 

minimum temperature occurs at the top. 

Figures 4.27 and 4.28 are plotted in one graph in Figure 4.29. In this fig- 
ure, one can see the sudden temperature drop and sudden temperature rise of 
the top portion of the core at the beginning of the shade and sun time respec- 
tively denoting the sensitivity of the core of the air to the boundary side wall 
temperature at different axial locations. 

Figure 4.30 shows the axial distribution of the local convective heat transfer 
coefficient vs. axial distance from the receiver aperture during the sun time. The 
solid line represents the h 2 at the beginning of the sun time and the dashed line 
denotes the h z at the end of the sun time. The following can be noticed from 
this figure: 
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1. ) At the beginning of the sun time minimum h z value is obtained at the 

middle of the core, however, at the end of the sun time the minimum value 

is obtained at the bottom of the core. 

2. ) Range of the h z values are much higher in the lower portion of the core. 

3. ) h z is decreasing during the sun time at all axial locations. 

Also, Figure 4.31 shows cyclic variation of the h for three different axial 
locations, namely, at z/L=.l , .5 and .95. Also, Figure 4.31 shows cyclic variation 
of the h for different axial locations. In the present study, It was experienced 
that, during the shade time, h z — ► oo due to the very small value of calculated 
boundary layer thickness. To avoid numerical difficulty, minimum boundary 
layer thickness has been set to a certain value which results in h as 70 Btu/hr ft 
°F which is constant during most shade period. 

Similarly, local Nusselt number ditribution is plotted in Figures 4.32 and 4.33. 
In Figure 4.32 the axial distribution of the Nu* vs. axial distance from the re- 
ceiver aperture is shown during the sun time. The solid line represents the Nu 2 
at the beginning of the sun time and the dashed line denotes the Nu z at the end 
of the sun time. The following can be noticed from this figure: 

1. ) It is shown that almost identical shape is maintained at the beginning of 

the sun time and at the end of the sun time. 

2. ) Nu z increases sharply near the mixing region indicating that more convec- 

tive heat transfer is occuring at the surface near the top, i.e., the outlet 

section of the receiver. 

Also, Figure 4.33 shows cyclic variation of the Nu for the three different 
axial locations mentioned in Figure 4.31. Due to the same limitation on the 
boundary layer thickness, constant limitting values are shown in this plot during 
most shade period. 
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In order to estimate the effects of the free convection heat transfer on 
the temperature of the receiver, the results with natural convection (NC) and 
without NC are compared in Figures 4.34 to 4.45. 

Figure 4.34 shows the temperature variation of the receiver components 
located at the inlet section of the working gas vs. time with natural convection 
included. In order to estimate the effect of natural convection on the receiver 
temperatures at the above locations, similar results to the Figure 4.34 is shown 
in Figure 4.35 with no natural convection included. Comparing Figures 4.34 
and 4.35, the following can be noticed: 

1. ) During the initial transient start-up cycles, temperatures of the receiver 

components have been reduced due to convection loss to air inside the 
cavity. 

2. ) The temperature of the working gas at this section remained nearly same 

as the gas temperature obtained from the case without natural convection. 

3. ) The maximum temperature of TES containment during the last cycle has 

dropped noticeably due to natural convection heat loss during sun time. 
Also, the average temperature of the TES containment at this section dur- 
ing the last cycle has dropped due to convection loss during sun time. 

4. ) During the last cycle, the temperature range of TES containment in this 

section has dropped about 50 °F due to the presence of the natural con- 
vection. 

Similarly, Figure 4.36 shows the temperature variation of the receiver com- 
ponents at the middle section of the working gas vs. time with natural convection 
involved. The effect of natural convection on this section of receiver temperatures 
can be observed by comparing this figure with Figure 4.37 which shows similar 
results with no natural convection included. Comparing Figures 4.36 and 4.37, 
the following can be noticed: 
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1. ) Again, during the initial transient start-up cycles, temperatures of the 

receiver components have been reduced due to convection loss to air inside 
the cavity. 

2. ) The maximum temperatures of the TES containment during the last cycle 

has dropped due to convection heat loss during the sun time. However, 
the minimum temperature of TES containment has increased due to heat 
supply from the air during the shade time. 

3. ) The average temperatures of the TES containment at this section during 

the last cycle has increased mainly due to heat supply during the shade 
time. (The heat supply during the shade time had more effect than the 
heat loss during the sun time.) While the average temperatures for the 
other receiver components remained nearly unchanged. 

4. ) The temperature ranges of these receiver components have been reduced 

about 27% due to natural convection. 

Similarly, Figure 4.38 shows the temperature variation of the receiver com- 
ponents at the outlet section of the working gas vs. time with natural convection 
included. Again, the effect of natural convection can be observed by comparing 
this figure with Figure 4.39 which shows similar results with no natural convec- 
tion included. Comparing Figures 4.38 and 4.39, the following can be noticed: 

1. ) Again, during the initial transient start-up cycles, temperatures of the 

receiver components have been reduced due to convection loss to air inside 
the cavity. 

2. ) Unlike the other sections, i.e., inlet and middle sections, the maximum 

temperature of the TES containment during the last cycle at the outlet 
section remained nearly unchanged. Also, the minimum temperature of 
TES containment has been increased noticeably due to natural convection. 
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3. ) The average temperature of the TES containment has increased mainly 

because of the heat supply from the air during the sun time. Due to 

relatively low maximum temperature of the TES containment, heat supply 

occurs during more than half of the sun time. 

4. ) The temperature ranges has been dropped to about 30% due to natural 

convection. 

Solid-liquid interface location vs. time is plotted in Figure 4.40 at three 
different axial nodes located at the inlet, middle and outlet section of the working 
gas tube with natural convection. Again, the effect of natural convection on the 
TES containment temperature can be observed by comparing this figure with 
Figure 4.41 which shows similar results with no natural convection included. 
From the Figures 4.40 and 4.41, the following can be noticed: 

1. ) Due to natural convection, more PCM is in liquid state in the middle and 

outlet section of the gas tube at the end of shade time. 

2. ) Due to natural convection, during the sun time, more PCM is melted at 

the fluid outlet section than the PCM at the middle section. 

Also, two dimensional motion of the interface vs. axial distance is shown in 
Figures 4.42 and 4.43 during sun period and shade period repectively. The effect 
of natural convection can be observed by comparing these figures with Figure 
4.44 and 4.45 which show similar results with no natural convection included. 
Comparing Figures 4.42 and 4.43 with Figures 4.44 and 4.45, the following can 
be noticed: 

1. ) At the inlet section, less PCM melts, however, more PCM melts at the 

outlet section of the gas tubing due to natural convection. 

2. ) The PCM solidifies more rapidly during the shade time than it melts during 

the sun time. 
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3.) At the end of shade time, most PCM is still in liquid state except the the 

fluid inlet section with natural convection included. 

Also comparison has been made between the result with natural convection 
(NC) analysis and without NC analysis and the temperatures of the important 
receiver components are tabulated in Table VI. 

The plots of figures 4.34 to 4.45 are su mm arized in Table VI to give more 
elucidatory comparison between the terrestrial and orbital results. Effect of 
thermal stratification has been detected from this enclosure natural convection 
analysis even though negative temperature gradient has occured during the shade 
time. As a result of the effect, the temperature of receiver components at the 
outlet section of gas tube has been increased rapidly and is maintaining higher 
average temperature through the cycle than the results obtained from the orbital 
test which has no natural convection effect (Table Via). 

The most noticeable outcome from this comparison is that the time aver- 
aged temperatures of the receiver components during terrestrial test are almost 
equal to those obtained from orbital test, instead of giving reduced temperature 
as in the previous preliminary results. From the results obtained, it is found that 
the existance of the air and resulting free convection does not have considerable 
effects on the receiver temperature with the receiver in vertical position with 
aperature pointing downwards. 

This is explained as follows by comparing with the preliminary results dis- 
cussed in previous section. In the preliminary study, the temperature of the core 
region inside a cavity was considered as ambient temperature or reservoir tem- 
perature which is constant all the time, and no variation in core temperature was 
assumed. As a result, substantial amount of the convective heat is lost contin- 
uously and this heat loss causes the overall decrease in in receiver temperature. 
On the contrary, the core temperature in the detailed study discussed in this 
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section, was calculated neither same as the ambient temperature nor constant 
all the time during a cycle. Also, with shade time included, the heat loss during 
the sun time is balanced with the heat supply from the air during the shade time 
due to sudden change of the wall temperature which is below the air temperature 
(see Figure 4.26). 

As discussed in Chapter III, the heated fluid near the wall moves upward 
and mixed together in the upper portion of the cavity and comes down through 
the core and then the fluid is entrained into the boundary layer again. Since 
the flow is circulating inside a cavity and has little interaction with the ambient 
temperature based on the assumption made previously for the aperture opening, 
incoming solar energy is conserved inside a receiver cavity without any heat loss 
to the ambient air outside the aperture. Accordingly no continuous convective 
heat loss is shown from terrestrial free convection results except during warm-up 
time of the air. 

It is worth noting that the range between the minimum temperature and 
the maximum temperature during a cycle is reduced due to natural convection, 
and one can see the thermal behavior of the air inside a cavity is very similar to 
that of energy storage material, i.e., cooling the receiver during sun time (storing 
the energy) and heating the receiver during shade time (liberating the energy). 
Therefore one should expect the receiver thermal performance would be worse in 
the orbit than the terrestrial results with free convection effects which predicts 
less temperature variation and lower maximum temperature of each component 
of the receiver which might ignore the material failure due to larger thermal 
expansion and/or higher maximum temperature than expected during the orbit 
operation. 


Table VI. Comparison of Lhc terrestrial TES receiver simulation 
results with the orbital TES receiver simulation results 



* Averaged values through the cycle 
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Figure 4.6 Working fluid temperature variation with time at dif 
ferent axial location, (with thermal energy storage) 
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Figure 4.7 Working fluid temperature variation with time at dif- 
ferent axial location, (without thermal energy storage) 
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Figure 4.9 The temperature variation of heat exchanger tubing 
with time at different axial location, (without thermal energy stor- 
age) 
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Figure 4.10 The temperature variation of the TES containment 
with time at different axial location, (with thermal energy storage) 
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Figure 4.11 The temperature variation of the inside receiver sur- 
face with time at different axial location, (without thermal energy 
storage) 
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Figure 4.12 Dimensionless interface location as function of time 
at different axial location. 
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Figure 4.13 The variation of liquid percentage of the thermal en- 
ergy storage material with time. 
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Figure 4.14 Two-dimensional interface motion during sun time for 
different times. 
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Figure 4.16 Energy transferee! in and out of the receiver as a func- 
tion of time. 
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Figure 4.18 Temperature distribution of the cavity surface using 
h = 6.86 W/m 2 • K. 




TURE 





ELAPSED TIME 


Figure 4.22 The core temperature variations of the air with time 
at z/L=QA. 
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Figure 4.27 Axial variation of the air temperature during the sun 
period, 7.825 < t(hr) < 8.792. 
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Figure 4.30 Axial distribution of local convective heat transfer 
coefficient at the beginning and the end of the sun time, t(hr) = 7.83 
and t{hr) — 8.79. 
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Figure 4.31 Cyclic variation of the convective heat transfer coef- 
ficient at three different axial locations, z/L= 0.1, 0.5 , and 0.95. 
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Figure 4.32 Axial distribution of local Nusselt number at the be- 
ginning and the end of the sun time, t(hr) = 7.83 and t(hr) = 8.79. 
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Figure 4.33 Cyclic variation of the Nusselt number at three dif- 
ferent axial locations, z/L— 0.1, 0.5 , and 0.95. 
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Figure 4.34 Temperature variation of the receiver components 
with time at the inlet section of the working gas. (with natural 
convection) 
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Figure 4.37 Temperature variation of the receiver components 
with time at the middle section of the working gas. (without natural 
convection) 
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Figure 4.39 Temperature variation of the receiver components 
with time at the outlet section of the working gas. (without nat- 
ural convection) 
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Figure 4.40 Dimensionless interface location as function of time 
at different axial location, (with natural convection) 
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Figure 4.41 Dimensionless interface location as function of time 
at different axial location, (without natural convection) 
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Figure 4.42 Two-dimensional interface motion during sun time for 
different times, (^th^nafural convection) 
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Figure 4.43 Two-dimensional interface motion during shade time 
for different times, (with natural convection) 
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Figure 4.44 Two-dimensional interface motion during sun time for 
different times, (TOthout natural convection) 
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CHAPTER V 

SUMMARY AND CONCLUSIONS 

fi-1 Summary 

The Space Station Freedom solar dynamic power generating system con- 
cept utilizes solar energy as the heat source in a closed Brayton cycle to develop 
auxiliary electrical power in space. In the present work, attention has been given 
to the solar receiver. The solar receiver is a combination of heat exchanger and 
thermal energy storage device which transfers a certain amount of heat to a 
Brayton cycle working gas during both sun and shade periods of an earth orbit. 
During a sun period, a parabolic collector has been used to focus solar radia- 
tion into the receiver through an aperture located at the collector focal point. 
The solar energy input is normally designed to be greater than that required for 
Brayton power system operation during a sun period. Now the excess energy 
can be stored within the receiver energy storage device as heat of fusion of phase 
change material (PCM) and then can be withdrawn during the following shade 
period to maintain power system operation. 

Heat transfer analysis of phase change material has been made to iden- 
tify the accurate location of the solid-liquid interface and this analysis is then 
applied to the the receiver thermal energy storage subsystem to evaluate the var- 
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ious parameters involving heat transfer through the PCM. Because of extreme 
complexity of the PCM thermal behavior, such as natural convection due to 
density difference between the two phases and void formation due to thermal 
expansion, present solution scheme is confined to heat conduction analysis with- 
out void or natural convection effects. However, this simplified thermal model 
makes the analysis tractable and yet retains the most important characteristics 
of the receiver. 

Due to time and cost constraints it is likely that most components of the 
space station solar dynamic power system will undergo extensive earth-based 
testing as opposed to orbital testing. This will expose them to gravity driven 
effects due to the presence of the gravity and the air. The presence of gravity 
and the air has been identified as the parameter that can significantly alter the 
performance of the receiver during a 1-g test operation of the system. 

Hereby, in the present work, another attention is given to the natural con- 
vection inside a receiver cavity to estimate its effects on the receiver performance 
during 1-g operation. Due to the lack of appropriate emperical formula for in- 
ternal natural convection problem which can be assumed to be suitable for the 
receiver cavity geometry, more rigorous method has been presented which in- 
cludes entire flow field inside a cavity and identify different flow characteristics 
within the domain. 

As a first step to the proposed receiver thermal analysis described above, 
overall energy balance involved in the solar receiver has been studied to under- 
stand conceptual design objective and scope in Chapter I. Also, the method of 
approximating amount of energy requirements from the other given energies was 
obtained based on lumped capacitance assumption. 

In Chapter II, the method of predicting thermal behavior of axisymmet- 
rical phase change material resides between two concentric cylinder, has been 
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investigated using non dimentionalized variables. Both implicit and explicit nu- 
merical scheme was tested using different number of radial nodes of 4, 6, 9, 12, 
15, 18 for the one dimentional ice formation problem which has been solved an- 
alytically. The location of the solid/liquid interface of the material as well as 
the temperature distribution of the phase change material has been calculated. 
From the results, implicit scheme has been identified unstable for this problem 
and explicit scheme gives stable results and the results are in good agreement 
with the analytical results. After estimating required CPU time and solution 
accuracy the explicit scheme with 9 radial nodes was selected and used in anal- 
izing receiver TES heat transfer. Instability of the implicit scheme has been 
presumed to result horn the increasing round off error due to large number of 
arithmatic operations required to invert the matrix for the increasing number of 
radial nodes. As an alternative, solving the matrix using iteration method which 
is not used in present study, may overcome the difficulty. Or using alternate di- 
rection implicit (ADI) method which makes the matrix tridiagonal is considered 
worth to try to improve the results of the simple implicit scheme. 

In Chapter III, the flow inside a closed cylindrial enclosure with side wall 
heated has been investigated using similar solution algorithm developed by Hess 
and Miller to analyze natural convection effects inside a receiver cavity during 
1-g operation. The flow domain was divided into three regions (Figure 3.1), 
namely, boundary layer, mixing, core region as suggested by previous investiga- 
tors. In each region, different method has been applied according to the flow 
characteristic. The velocity and temperature distribution are calculated in the 
boundary layer region using Karman-Pohlhausen method and the temperature 
of the mixing region is obtained by performing the energy balance over the con- 
trol volume of mixing region. Also one dimensional velocity and temperature 
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has been calculated in the core region. The three regions are coupled with each 
other, therefore a solution of one region becomes the boundary condition of the 
other region. 

Numerical data obtained for water inside a circular cylinder subject to two 
different boundary temperature at the top perimeter of the cylinder are compared 
with Miller’s (1977) numerical solution which uses full Navior-Stokes equation 
and Hess’s integral results. Velocity distribution inside the boundary layer and 
temperature distribution of the core region are in good agreement with Miller’s 
results, and due to new temperature profile modefied from the parabolic profile 
suggested by Hess, excellant agreement were shown in predicting temperature 
distribution in mixing region which has been considered as important region 
connecting boundary layer and core region. 

The numerical method of Chapter II and III has been applied to the Bray- 
ton cycle solar dynamic receiver thermal analysis. And two subroutines devel- 
oped in Chapter II and III, for TES analysis and natural convection analysis 
respectively, have been incorporated into a computer program written for anal- 
izing solar receiver thermal performance which does not have capability of those 
analyses. In Chapter IV, the results are compared with the results obtained 
without TES analysis and/or natural convection analysis. 

5.2 Conclusions 

5.2.1 Effect of TES on Receiver Thermal Performance 

It has been identified that the temperature variation of the working gas 
at the outlet section during an orbit cycle can be reduced more than 50 % by 
adding phase change material to the place between the inside cavity surface and 
working gas tube so that more constant power output of the entire system can 
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be obtained which is one of the most important objective of the solar dynamic 
power system. However, higher maximum surface temperature resulting from 
thermal energy storage, has been observed which might cause the damage or 
structural failure of the TES containment material. 

During the sun time, incoming solar energy has been used to increase 
the surface temperature due to larger thermal resistance of the PCM while the 
energy is transfered to the working gas at much faster rate when no TES is 
present. As a result, different effect has been identified for different receiver 
components according to relative location of the receiver component to PCM of 
thermal energy storage. For working gas and gas tubing, temperature fluctuation 
during a cycle has been significantly reduced (60%) with almost same average 
temperature. For cavity surface, on the other hand, almost same temperature 
fluctuation was obtained with noticeably increased average temperature. 

Since only 60 % of the PCM has been melted during the cycle and most 
PCM solidifies completely at the end of the shade time, it is desirable to increase 
the amount of the liquid PCM during the sun time in order to store more energy. 
If mass flow rate of the working gas is reduced during sun time then more heat 
would be stored in PCM due to reduced heat transfer to working gas. 

5.2.2 Effect of Natural Convection on Receiver Performance 

On the natural convection effects, nearly the same temperatures for the 
overall receiver component, including gas outlet temperature, has been obtained 
in terrestrial results as seen in orbital results. The effect of thermal stratification 
has been identified which forms the temperature gradient in the axial direction. 
As a result, local temperature increase at the outlet section and temperature 
drop at the inlet section of the gas tube was obtained. 

The most noticeable outcome from this comparison is that the time aver- 


154 


aged temperatures of the receiver components during terrestrial test are almost 
equal to those obtained from orbital test, instead of giving reduced temperature 
as in the previous preliminary results. From the results obtained, it is found that 
the existance of the air and resulting free convection does not have considerable 
effects on the receiver temperature with the receiver in vertical position with 
aperature pointing downwards. 

This is explained as follows by comparing with the preliminary results dis- 
cussed in previous section. In the preliminary study, the temperature of the 
core region inside a cavity was considered as ambient temperature or reservoir 
temperature which is constant all the time, and no variation in core tempera- 
ture was assumed. As a result, substantial amount of the convective heat is lost 
continuously and this heat loss causes the overall decrease in receiver temper- 
ature. On the contrary, the core temperature in the detailed study discussed in 
this section, was calculated neither same as the ambient temperature nor con- 
stant all the time during a cycle. Also, with shade time included, the heat loss 
during the sun time is balanced with the heat supply from the air during the 
shade time due to sudden change of the wall temperature which is below the air 
temperature. 

It is worth noting that the range between the minimum temperature and 
the maximum temperature during a cycle is reduced due to natural convection, 
and one can see the thermal behavior of the air inside a cavity is very similar to 
that of energy storage material, i.e., cooling the receiver during sun time (storing 
the energy) and heating the receiver during shade time (liberating the energy). 
Therefore one should expect the receiver thermal performance would be worse in 
the orbit than the terrestrial results with free convection effects which predicts 
less temperature variation and lower maximum temperature of each component 
of the receiver which might ignore the material failure due to larger thermal 
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expansion and/or higher maximum temperature than expected during the orbit 
operation. 


5.3 Future Work 


Although the present model incorporates complicated aspects of the heat 
transfer and fluid mechanics in a consistant manner, the following areas can be 
further developed for an improved and more comprehensive model: 

1. ) The formation of the void in the phase change material (PCM) due to 

thermal contraction should be included in the analysis to determine more 
accurate heat resistance of the PCM . 

2. ) The convection heat transfer inside the liquid PCM should be included 

by allowing density variation between the solid and liquid PCM to obtain 
better approximation on the heat transfer phenomena that takes place in 
actuality. 

3. ) A reliable and stable linearizing procedures should be utilized to overcome 

the numerical difficulties experienced in solving coupled non-linear bound- 
ary layer PDE to allow larger time step size in implicit scheme. 

4. ) The properties variation of the receiver material due to high temperature 

range (300 °F maximum) during the orbit cycle, should be included to 
estimate more realistic heat transfer. 

Despite some of the limitations of the present heat transfer and flow model, 
it is believed that the essential features of the phase change phenomena of the 
thermal energy storage system and flow fields inside a receiver cavity have been 
illustrated. The present work should be viewed as a first logical step towards a 
more detailed analysis of the phase change and enclosure free convection problem. 
In this regard, the future developments pointed above would provide a better 


understanding for these problem. 
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APPENDIX 

FOURIER (von NEUMANN) STABILITY ANALYSIS 
FOR 2-D AXISYMMETRIC HEAT EQUATION 

A finite-difference approximation to a partial differential equation may be 
consistent but the solution will not necessarily converge to the solution of the 
PDE. According to the Lax Equivalence theorem [see Anderson, 1984] a stable 
numerical method must also be used. Purpose of the analysis in this section is 
to derive a conditional atability requirement on the time step and the spatial 
mesh spacing for the heat equation which was treated in Chapter II via Fourier 
or von Neumann stability analysis. 

Consider the two-dimensional axisymmetric heat equation 
(see Equation 2.2): 


du p ( d 7 u ^ 1 du ^ 


dt 


^ dr 2 r dr "** dz 2 ) 


+ 


( 1 ) 


The simple explicit finite-difference approximation for this problem is: 


_ -,n 

u j,fc u i,fc 

At 
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( 2 ) 
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Or 


where 


•JS 1 ~ U l,k + r l ( u j+l,k - 2u j,k + u j- l,k) + r 2 ( u j+l,k - u j-l,k) 
+ T * ( U j,k+1 ~ 2u j,k + u j,k-l) 

rj = 0(Ai/Ar 2 ) 
r 2 — 0(At/2jAr 2 ) 
r 3 = 0(At/Az 2 ) 


( 3 ) 


( 4 ) 


In this case a Fourier component of the form is assumed: 


„n _ -at ik T r ik,z 
u j,k ~ e e e 


( 5 ) 


where k r and k z are the wave numbers for r and z coordinates respectively. 
If Equation (5) is substituted into Equation (3) we obtain 

e a(t+At) e ik T r e ik,z _ e at £ ik r r £ ik t z + ^ + ^at £ ik r (r+ Ar) £ ik M z 

+ (rj - r 2 )e at e ikr(r ~ Ar) e ik ' z 


-2 (n +r 3 )e al e ik ' r e ik ' z 
+ r 3 e at e ik,(z+Az) 


+ r 3 e at e ik rr e ik t {z-Az) 


if we divide by e at e xkrT e xktZ and utilize the relation 


( 6 ) 


a e ip + e xp e x P - e 

cos p = and sin p = 


the above expression in Equation (6) becomes 
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e aAl = 1 + (n + T 2 )e ik ' Ar + (n - r 2 )e-^ Ar - 2(r 3 + r 3 ) 

+ r*(e <4 ' A * + e- i ** A *) 

= 1 + r 1 (e < *' Ar + c _tt ' Ar ) + r 2 (e ifc ' Ar - e" ifc ' Ar ) 

+ r 3 (e rt - A * + e~ ik ‘* z ) - 2 pj - 2r 3 

e ik T Ar _|_ g— *fc,Ar ^ ^ gifc, Ar _ £ —ik r Ai* ' 


/ e **rAr , -*fc r Ar \ / Ar _ -»fc r Ar \ 

= l + 2r,( j 1 }+^( 2 ) 

/ Az i --tfc, Az \ 

(—2 >) 


+ 2r s 


= 1 + 2r 1 (cos/?i — 1) + 2r 3 (cos/? 2 — 1) + 2 r 2 isin/ 3 x 
where fix = fc r Ar and / 3 2 = k z Az. Employing the trigonometric identity 


. , /? 1 — cos /? 

in z — = 


sin 


the amplification factor is 


G = e aAt = 1 — 4 7*! sin 2 — ■ — 4 r 3 sin 2 — 4. 2r 2 i sin fix 


Thus for stability; 


where 


\G \<1 


( 7 ) 


( 8 ) 


|C7| = |l — 4 t*i sin 2 — — 4 r 3 sin 2 — + 2r 2 i sin 0x \ 

= y^l— 4rjsin 2 — 4r 3 sin 2 ^ | + (2r 2 1 sin /?! ) 3 


( 9 ) 


After rearranging the Equation ( 9 ), the stability requirement Equation (8) is 
then 

1 1 £ 2 *4 ^ 1 1 1 2 *4 ^2 q *2 o • 2 ^^2 

1 + 16 rj sm — + lorj sm — — Sri sm — — 8 r 3 sm — 
z z z z 

32 rir 3 sin 2 ^ sin 2 ^ + 4 r| sin 2 fix <1 


(10) 
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Further reduction yields 

4r 2 -f- 4r| — 2 rj — 2t*3 + 8 rir 3 + t\ < 0 


Substituting , r 2 , and r 3 from Equation (4) gives 

4/? 2 Af 2 4,9* Ai* J?A< ftAf <3 2 A < 2 £ 2 At 2 

Ar* + A 2 * Ar 2 Az 2 + Ar 2 Az 2 + 4Ar 2 “ 

which yields the following constraint on the size of the time step relative to the 
size of the mesh spacing: 


8 Ar 2 Az 2 (Ar 2 + Az 2 ) 

“ /? 16Ar 4 + 16 Az 4 + 32Ar 2 Az 2 + Ar 2 Az 4 


( 11 ) 


The Equation ( 11 ) , which is the same equation as Equation 2.3 in Chap- 
ter II, has been used for explicit method of phase change problem in the chapter. 


